**3. Investigation of drug resistant mechanism by MD simulation**

Genetics study of influenza H5N1 virus isolated from patients who died despite being given oseltamivir showed that mutations H274Y or N294S confer high-level resistance to oseltamivir (De Jong et al., 2005 & Le et al., 2005)). There is even emerging evidence that these drug-resistant mutants pose the same risk with H1N1pdm, as shown by reported cases of the H274Y mutation of H1N1pdm (Guo et al, 2009). The rapid emergence of oseltamivir resistance in avian flu has already motivated numerous studies, both experimental and theoretical, to uncover the mechanisms of how point mutations in neuraminidase alter drug binding. Despite initial inroads, the current understanding of drug resistance remains incomplete and some conclusions are conflicting. For example, in one study, it was reported that the H274Y mutation disrupts the E276-R224 salt bridges, but in a separate study the same salt bridges were observed to be stable. We characterize the drug-protein interactions of oseltamivir bound forms of wild type avian H5N1 and swine H1N1pdm neuraminidases and how their mutations, H274Y and N294S, rupture these interactions to confer drug resistance through molecular modelling and atomistic MD simulation.

### **3.1. Computational details**

166 Bioinformatics

Rank / Comp. ID

NSC

standard deviations.

parallel automatic pipeline.

HM mean energy

(kcal/mol)

Predicted Ki (

7 Oseltamivir -8.69 0.429 -8.57 0.35

AM mean energy

(kcal/mol)

**Table 1.** Docking results ranked by harmonic mean of binding energy. The first column is the final rank and also the compound ID. The predicted Ki calculated according to the harmonic mean binding free energies are also shown, as well as the arithmetic mean binding free energies and their corresponding

The results also confirm the observed effectiveness of oseltamivir against both the swine H1N1pdm and H5N1 flu. Detailed analysis on the hydrogen bond networks between top binding candidates with the swine H1N1pdm neuraminidase protein reveals that the mutations H274Y and N294S do not have any direct interactions with these compounds, and thus suggest the possibility of using the present top hit compounds for further computational and experimental studies to design new antiviral drugs against swine H1N1pdm flu virus and its variants (Hung et al., 2009). Furthermore, a more complete virtual screening using the full NCI diversity set (NCIDS) or larger sets from the ZINC database should be done to identify drug candidates that may have been missed in this study. The NCIDS has over 2000 compounds thus the full screening should be done with

**3. Investigation of drug resistant mechanism by MD simulation** 

Genetics study of influenza H5N1 virus isolated from patients who died despite being given oseltamivir showed that mutations H274Y or N294S confer high-level resistance to oseltamivir (De Jong et al., 2005 & Le et al., 2005)). There is even emerging evidence that these drug-resistant mutants pose the same risk with H1N1pdm, as shown by reported

Standard deviation

(kcal/mol)

Structure

μM)

The coordinates for the H5N1 neuraminidase bound with oseltamivir was taken from a monomer of the Protein Data Bank (PDB) structure 2HU4 (tetramer), while those of mutants H274Y and N294S were taken from structures 3CL0 (monomer) and 3CL2 (monomer) respectively. Even though the biological form of neuraminidase is tetrameric, its monomer contains a functionally complete active site and yields reasonable results in a prior study using MD simulations. The position for oseltamivir bound to H1N1pdm was adopted from its corresponding location in H5N1, as the two proteins' binding pockets differ only by residue 347 (which is Y in H5N1 and N in H1N1pdm), located on a loop at the periphery of the active site. Oseltamivir-mutant complexes of H1N1pdm were built by mutating H274Y and N294S of the H1N1pdm wild type model. In total, 6 systems were modelled and simulated for oseltamivir bound H5N1, and H1N1pdm wild type and H274Y and N294S mutants. Crystallographically resolved water molecules and a structurally relevant calcium ion near the native binding site for SA were retained and modelled in all simulated systems. The protein complexes were then solvated in a TIP3P water box and ionized by NaCl (0.152M) to mimic physiological conditions.

All simulations were performed using NAMD 2.7 and the CHARMM31 force field with CMAP correction. The ionized systems were minimized for 10,000 integration steps and equilibrated for 20 ns with 1 fs time steps. Following this, a 20 ns unconstrained equilibration production run was performed for subsequent trajectory analysis, with frames stored after each picosecond (every 1000 time steps). Constant temperature (T = 300 K) was enforced using Langevin dynamics with a damping coefficient of 1 ps−1. Constant pressure (p = 1 atm) was enforced using the Nosé-Hoover Langevin piston method. Van der Waals interaction cut-off distances were set at 12 Å (smooth switching function beginning at 10 Å) and long-range electrostatic forces were computed using the particle-mesh Ewald (PME) summation method.

Analysis included the calculation of an averaged electrostatic potential field over all frames of the trajectory using RMSD-aligned structures. Maps of the electrostatic potential field were calculated on a three-dimensional lattice. The long-range contributions to the electrostatics were approximated using the multilevel summation method (MSM), which uses nested interpolation of the smoothed pairwise interaction potential, with computational work that scales linearly with the size of the system. The calculation was performed using the molecular visualization program VMD that provides a graphics processing units (GPU) accelerated version of MSM to produce the electrostatic potential map. The GPU acceleration of MSM ( Hardy et al., 2009) provided a significant speedup over conventional electrostatic summation methods such as the Adaptive Poisson Boltzman Solver (APBS), achieving a benchmark processing time of 0.2s per frame versus 180 seconds per frame (on a conventional CPU) using APBS69 for a 35,000 atom system, offering a speedup factor of about 900. The use of GPU acceleration enabled averaging the electrostatic potential field over all frames of the simulation trajectories. Hydrogen bond analysis utilized a distance and angle cut-offs of 3.5 Å and 60 degrees, respectively.

Incorporating Molecular Dynamics Simulations

into Rational Drug Design: A Case Study on Influenza a Neuraminidases 169

strains significantly alters the drug-protein stability in regard to the hydrogen bond network involved. H274Y mutation induced disruption of the stable hydrophobic packing of

**Figure 3.** Network and occupancy of hydrogen bonds stabilizing oseltamivir in the SA binding pocket of wild type and drug-resistant mutant avian H5N1 neuraminidases, in simEQ1, simEQ3, and simEQ5.(A) shows histograms of the percent of hydrogen-bond occupancies for interactions between oseltamivir and residues E119, D151, R152, R292, Y347, and R371 across each simulation run. (B) through (D) are schematic views depicting the orientation of protein sidechains which form protein-

drug hydrogen bonds.

oseltamivir's pentyl group in both H5N1 and H1N1pdm neuraminidases.
