**2.2. Ligands**

The high percentage sequence identity (91.47 %) of swine H1N1pdm compared to avian H5N1 neuraminidase, and similarities in 3D structures (active-site 150-loop and SA binding site residues) explain why oseltamivir, which has been developed for the H5N1 virus, is also effective against the swine flu virus H1N1pdm. These facts also suggest that top-binding ligands of avian H5N1 are likely to have high affinity to the swine flu H1N1pdm neuraminidase. Therefore, the 27 top hits for H5N1 neuraminidase (Cheng et al., 2008) and six known drugs for avian H5N1 namely: SA, N-acetyl neuraminic acid, aka, NANA or Neu5Ac), DANA (2,3-didehydro-2-deoxy-N-acetyl neuraminic acid), oseltamivir, zanamivir, peramivir (clinical trial phase 3 drug candidate) and SK (shikimic acid) were selected for our study. Ligand structures were obtained from the NCI and optimized by the MOPAC (Stewart, 2007). Our selection of ligands for docking is based on the high sequence identity and observation that mutations mainly appear on the H1N1 neuraminidase surface not its active site (Le et al., 2009). However, further screening for all the entries of the NCI is suggested to avoid false positives and false negatives.

#### **2.3. Molecular docking**

162 Bioinformatics

in cyan.

**2.2. Ligands** 

representative clusters were obtained, which account for 96.2% of the configuration space from

**Figure 2.** 13 representative ensembles resulting from clustering analysis, accounting for 96.2% of the configuration space, are ordered by the corresponding simulation time. The four largest cluster ensembles account for over 66% of the configuration space from the 20 ns MD simulation. The most dominant cluster is colored in orange, the second one in silver, the third one in blue and the fourth one

The high percentage sequence identity (91.47 %) of swine H1N1pdm compared to avian H5N1 neuraminidase, and similarities in 3D structures (active-site 150-loop and SA binding site residues) explain why oseltamivir, which has been developed for the H5N1 virus, is also effective against the swine flu virus H1N1pdm. These facts also suggest that top-binding ligands of avian H5N1 are likely to have high affinity to the swine flu H1N1pdm neuraminidase. Therefore, the 27 top hits for H5N1 neuraminidase (Cheng et al., 2008) and six known drugs for avian H5N1 namely: SA, N-acetyl neuraminic acid, aka, NANA or Neu5Ac), DANA (2,3-didehydro-2-deoxy-N-acetyl neuraminic acid), oseltamivir, zanamivir, peramivir (clinical trial phase 3 drug candidate) and SK (shikimic acid) were selected for our study. Ligand structures were obtained from the NCI and optimized by the MOPAC (Stewart, 2007). Our selection of ligands for docking is based on the high sequence identity and observation that mutations mainly appear on the H1N1 neuraminidase surface not its active site (Le et al., 2009). However, further screening for all the entries of the NCI is

suggested to avoid false positives and false negatives.

the 20 ns MD trajectories. These 13 clusters were used in the docking experiments.

In the docking experiments, the 13 most representative configurations were used as receptors. AutoDockTools 1.5.2 was used to add polar hydrogens, assign Gasteiger charges11 and create grid binding boxes. The volume of each grid box was 72 x 72 x 72, with the default 0.375 Å spacing. The binding box was positioned to encompass all three possible binding sites, namely the SA, 150 and 430 cavities (Cheng et al., 2008). AutoGrid version 4.2.1 was used to calculate the binding affinities using the following atom types: A (aromatic carbon), C, N, NA (hydrogen bond accepting N), OA (hydrogen bond accepting O), P, S, SA (hydrogen bond accepting S), Cl, HD (polar hydrogen) and e (electrostatics). AutoDockTools version 1.5.2 was also used to merge nonpolar hydrogens, add Gasteiger charges and visually set up rotatable bonds for each ligand via AutoTors. For screening larger ligand libraries, Autodock Vina (Trott et al., 2010), which run more efficiently than this version, is highly suggested. The Lamarckian genetic algorithm was used to do the docking experiments using AutoDock 4.2.1.16 Docking parameters were chosen to reproduce structures of 13 corresponding oseltamivir–neuraminidase complexes in the MD simulation. Other parameters are as follows: trials of 100 dockings, population size of 200, random starting position and conformation, translation step range of 2.0 Å, rotation step range of 50 degrees, maximum number of generations of 27000, elitism of 1, mutation rate of 2%, crossover rate of 80%, local search rate of 6%, 8 million energy evaluations, unbound model was "same as bound", and docked conformations were clustered with the tolerance of 2.0 Å RMSD. Docking results were sorted by the lowest binding energy of the most populated cluster in cases of convergence. In the case of no dominant cluster, docking results were visually analyzed using VMD (Humphrey et al., 1996) to choose the best binding pose. Hydrogen bond analysis utilized a distance and angle cut-off of 3.5 Å and 45 degrees, respectively.

The appearance frequency (AF) of important hydrogen bonds is calculated as follows:

$$AF = \frac{\sum\_{i=1}^{13} n\_i h\_i}{\sum\_{i=1}^{13} n\_i} \tag{1}$$

Where i is the index number of each ensemble; ni the size of each ensemble; and hi = 1 if hydrogen bond exists, and 0 otherwise. After the dockings, statistical calculations were performed to obtain the final binding energies for compound ranking, using the arithmetic mean (AM) and harmonic mean (HM) binding energies as defined in the previous study. The arithmetic means were calculated directly from the binding energies

$$AM = \frac{\sum\_{i=1}^{13} n\_i E\_i}{\sum\_{i=1}^{13} n\_i} \tag{2}$$

Where Ei is the binding energy of each ensemble with the standard deviation

$$SD = \sqrt{\frac{\sum\_{i=1}^{13} n\_i \left(E\_i - AM\right)^2}{\sum\_{i=1}^{13} n\_i}}\tag{3}$$

The harmonic means were calculated by first converting the binding energies into inhibition constant Ki

$$K\_i = e^{\frac{1000E\_i}{RT}}\tag{4}$$

Rank / Comp. ID

NSC

HM mean energy

(kcal/mol)

Predicted Ki (

3 109836 -9.74 0.072 -9.50 0.47

4 350191 -9.56 0.099 -8.77 0.68

5 164640 -8.95 0.274 -8.78 0.42

6 5069 -8.85 0.326 -8.23 1.11

AM mean energy

(kcal/mol)

Standard deviation

(kcal/mol)

μM)

Incorporating Molecular Dynamics Simulations

Structure

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

Where R is the Boltzmann constant and T is the temperature (298.15K). The harmonic means were calculated using below equation then converted back to HM binding energies.

$$\overline{K\_i} = \frac{\sum\_{i=1}^{13} n\_i}{\sum\_{i=1}^{13} \frac{n\_i}{K\_i(i)}} \tag{5}$$

#### **2.4. Results**

The obtained docking scores reveal that 6 compounds, specifically NSC211332, NSC141562, NSC109836, NSC350191, NSC1644640, and NSC5069, have higher binding affinity to the H1N1pdm neuraminidase protein than oseltamivir does.


Incorporating Molecular Dynamics Simulations into Rational Drug Design: A Case Study on Influenza a Neuraminidases 165


164 Bioinformatics

constant Ki

**2.4. Results** 

Rank / Comp. ID

NSC

HM mean energy

(kcal/mol)

<sup>13</sup> <sup>2</sup>

*n E AM*

*i i*

*n*

(3)

*RT K e <sup>i</sup>* (4)

(5)

Structure

13 1

1000 *<sup>i</sup> E*

1

*SD*

*i i i*

The harmonic means were calculated by first converting the binding energies into inhibition

Where R is the Boltzmann constant and T is the temperature (298.15K). The harmonic means

13

1 13

The obtained docking scores reveal that 6 compounds, specifically NSC211332, NSC141562, NSC109836, NSC350191, NSC1644640, and NSC5069, have higher binding affinity to the

*i i*

*n*

*i i i*

*n K i*

<sup>1</sup> ( )

were calculated using below equation then converted back to HM binding energies.

*i*

AM mean energy

(kcal/mol)

Standard deviation

(kcal/mol)

*K*

H1N1pdm neuraminidase protein than oseltamivir does.

Predicted Ki (

1 211332 -12.05 0.001 -11.73 0.61

2 141562 -10.14 0.037 -9.87 0.63

μM)

Incorporating Molecular Dynamics Simulations

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

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

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

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)

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

simulation.

**3.1. Computational details** 

(0.152M) to mimic physiological conditions.

summation method.

**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 standard deviations.

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 parallel automatic pipeline.
