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

1. Introduction

Nowadays, information about targeted diseases needed for structure-based drug design can be accessed easily. There are more than 80,000 coordinate entries of macromolecules revealed by X-ray, NMR, and Electron Microscopy techniques available in the Brookhaven Protein Data Bank. Many new target proteins whose 3D structures have not been solved experimentally can also be easily predicted by homology modelling with adequate reliability. Similarly, thousands of compounds including their drug-like properties are provided in free chemical databases such as the NCI Diversity set, ZINC, and the Drug Bank. Therefore, it is crucial to applied suitable tools to take advantages of these promising resources for rational drug design. Molecular dynamics (MD) simulation has proven itself as a complement to lab experiments to fill in the gaps in our knowledge between 3D structures of biological targets and their potential inhibitors (Rognan et al., 1998 & Jacob et al., 2011). MD can be simply understood as a simulation of physical movements of particles like atoms in atomistic simulation or a group of atoms in coarse-grained simulation. These particles are allowed to move and interact with each others. Their motions are described by Newton's equations, in which forces between the particles and potential energy are obtained from empirical data or quantum mechanical calculations. There are quite several MD simulation packages that are free for academics such as NAMD (Phillips et al., 2005), AMBER (Case et al., 2005), GROSMOS (Christen et al., 2005) and CHARMM (Brooks et al., 1983). The method was introduced by Alder and Wainwright in the late 1950's (Alder and Wainwright, 1957) for the study of the behaviour of simple liquids and the first protein simulation, bovine pancreatic trypsin inhibitor was done successfully twenty years later (McCammon et al., 1977). After that, MD simulation has been used wildly for investigating, interpreting and discovering the dynamics of biological macromolecules, not only fast internal motions but also slow conformational changes and even whole folding processes of some small proteins. MD simulation has become important complement to experimental procedures in structure-function determination of targeted diseases and rational design of their inhibitors. In the lead discovery stage, thousands of compounds from free chemical databases are docked against the rigid structures of receptors. Such screening processes are fast but are associated with a high number of false positives and false negatives. Flexible docking has improved the outcome of this process, but to a very limited degree. A recent method that significantly increases the accuracy of the lead identification process is called ensemble-based docking. This method combines the docking algorithm (Morris et al., 1998) with the dynamic structure of proteins taken from MD simulations as well as statistical analysis of binding energy (Amaro et al., 2008). Another important aspect of rational drug design is to understand all drug-protein interactions at an atomistic level and how these interactions are ruptured by certain mutations, leading to drug resistance. Thorough understanding of drug resistant mechanisms is crucial for the design of new agents effective against current drug-resistant strains. MD simulation is also an important tool in identifying drug binding pathway. However, this could be one of the most challenging tasks in computer-aided drug design due to the limitation of MD simulation time scale. One of the solutions is to combine electrostatic surface potential analysis with steered MD (SMD) simulation (Isralewitz et al., 2001). Forces can be applied to drug molecules in order to pull it along possible pathways predicted from electrostatic surface potential. The findings of drug binding pathway are highly important for investigation on non-active site mutations which possibly prevent drugs from entering binding site by rupture the structure of binding pathway as well as for the rational design of new drugs with good binding kinetics.

In this chapter we highlight the successes in using MD simulations to speed up the rational design of antiviral drugs. Due to their ability to mutate quickly, the viruses can easily survive people’s immune systems and resist current drugs. Currently, there are two pathogenic viruses, the 2003 avian flu (H5N1) virus( Ford Stephen et al., 2006) and the 2009 swine flu (H1N1pdm) virus (Butler et al., 2009), which share the same subtype 1 (N1) of neuraminidase. Neuraminidase, a glycoprotein component on a flu virus’ surface that cleaves the alpha-ketosidic linkage of sialic acid (SA) located on human cells to release new virions that then infect nearby healthy cells, is the most important target for antiviral drug design. Oseltamivir functions as a neuraminidase inhibitor that binds competitively to the Sialic Acid (SA) binding site (Laver et al., 2006). Unfortunately, oseltamivir resistance has been a big concern to public health since 2005 and thus continues effort to gain insight of drug resistant mechanism and designing new effective medicine become highly significant.

2. Ensemble-based docking

2.1. Receptor

2.1.1. MD simulation of ionized oseltamivir-neuraminidase complex

The molecular model of H1N1pdm used in MD simulation was built from the crystal structure of avian H5N1 and the amino acid sequence of swine flu H1N1pdm from GenBank Locus ID CY041156. All simulations were performed using NAMD 2.720 and the CHARMM31 force field with the CMAP correction. The ionized systems were minimized for 10,000 integration steps and equilibrated for 20 ns with a 1 fs time step. 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. The trajectories from our 20ns production run of EQ simulation of H1N1pdm neuraminidase with oseltamivir bound, were clustered as shown in Figure 1, and then used for ensemble based docking.

Figure 1.

Schematic representations of drug-bound simulation systems. Shown here as a representative example of simulated systems is H1N1pdm bound to oselvamivir. In A), the simulation system is shown in the full explicit solvation box with oseltamivir and the active-site calcium ion labelled. In B), oseltamivir is shown buried in the SA binding pocket of H1N1pdm rendered in surface view

2.1.2. RMSD clustering to extract receptor ensembles from an all-atom MD simulation

The holo system with oseltamivir removed was used for RMSD clustering and the docking experiments. Clustering analyses were performed on 20 ns MD trajectories using the g_cluster tool in the Gromacs package. In brief, snapshots at every 10 ps over the 20 ns simulation were recorded. 2000 resulting structures were superimposed, using all Cα atoms to remove possible rotational and translational movements of the whole system. We have visually verified that the binding-site residues of avian H5N1 neuraminidase, which cover the active site and are responsible for interaction with putative inhibitors, can also be used for the swine flu H1N1pdm neuraminidase. The four most populated clusters are shown in Figure 2. The RMSD clustering analysis was performed on this subset (117-119, 133-138, 146-152, 156, 179, 180, 196-200, 223-228, 243-247, 277, 278, 293, 295, 344-347, 368, 401, 402, and 426-441) using all-atoms (including side chains and hydrogen atoms) with the cut-off of 1.5 Å. A total of 13 representative clusters were obtained, which account for 96.2% of the configuration space from the 20 ns MD trajectories. These 13 clusters were used in the docking experiments.

Figure 2.

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 in cyan.

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

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{{\displaystyle \sum _{i=\text{1}}^{\text{13}}{n}_i{h}_i}}{{\displaystyle \sum _{i=\text{1}}^{\text{13}}{n}_i}}$$E1

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{{\displaystyle \sum _{i=\text{1}}^{\text{13}}{n}_i{E}_i}}{{\displaystyle \sum _{i=\text{1}}^{\text{13}}{n}_i}}$$E2

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

$$SD=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^{13}{n}_i{\left(E_i-AM\right)}^2}}{{\displaystyle {\sum }_{i=1}^{13}{n}_i}}}$$E3

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

$$K_i=e^{\frac{1000E_i}{RT}}$$E4

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{{\displaystyle \sum _{i=1}^{13}{n}_i}}{{\displaystyle \sum _{i=1}^{13}\frac{n_i}{K_i(i)}}}$E5

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.

Rank / Comp. ID	NSC	HM mean energy (kcal/mol)	Predicted Ki (μM)	AM mean energy (kcal/mol)	Standard deviation (kcal/mol)	Structure
1	211332	-12.05	0.001	-11.73	0.61
2	141562	-10.14	0.037	-9.87	0.63
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
7	Oseltamivir	-8.69	0.429	-8.57	0.35

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.

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

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.

3.2. Results and discussion

Hydrogen bonds which form the primary interactions between oseltamivir and H5N1/ H1N1pdm neuraminidases were shown in figure 4 and figure 5 accordingly. While the SA binding sites of H5N1 and H1N1pdm appear to differ mainly in the sequence of loop residue 347, it is not well understood whether antiviral drugs bind to each protein in the same manner, or if the drug resistant H274Y and N294S mutations disrupt critical hydrogen bonds. To address this question, the hydrogen bonds which form between oseltamivir and the residues lining the SA binding pockets of H5N1 and H1N1pdm were calculated for all simulation trajectories. R292 and R371 were observed to hydrogen bond with oseltamivir’s carboxyl moiety, and E119 and D151 with oseltamivir’s amino group (NH₃⁺). The H274Y mutation, however, appeared to disrupt the hydrogen bonding of oseltamivir’s acetyl group with R152, an interaction which was seen in the wild type and N294S systems for both simEQ1 and simEQ2. Prior analyses of crystallographic data alone suggested that Y347 forms a stable hydrogen bond with oseltamivir’s carboxyl group and is the source of oseltamivir-resistance in the N294S mutant.

Our MD simulations however, reflect statistics collected from long timescale simulations which produce a dynamic picture of molecular interactions in greater detail and resolution than can be seen from a static crystal structure. In our simulations oseltamivir’s carboxyl group primarily forms hydrogen bonds with R292 and R371, having little involvement with Y347. In fact, residue Y347 undergoes rotation to interact strongly with residue W295. Therefore, the speculation from previous studies, that the N294S mutation in the case of H5N1 actually destabilizes the hydrogen bonding between oseltamivir and Y347 to induce drug resistance, is not supported in our simulations.

The notable difference between H5N1 and H1N1pdm neuraminidases is the replacement of Y347 by N347 at the drug binding pocket. No conserved drug-protein hydrogen bond was observed for N347 in any of the three H1N1pdm simulations. Given the transient nature of even the N294S mutant induced hydrogen bond involving residue 347 in the case of H5N1, and the lack of interaction with residue 347 in any of the other simulations, it is highly unlikely that the single residue change (Y347 to N347), between the H5N1 and H1N1pdm 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 oseltamivir’s pentyl group in both H5N1 and H1N1pdm neuraminidases.

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.

Figure 4.

Network and occupancy of hydrogen bonds stabilizing oseltamivir in the SA binding pocket of wild type and drug-resistant mutant avian H1N1pdm neuraminidases, in simEQ2, simEQ4, and simEQ6. (A) are histograms of the percent of hydrogen-bond occupancies for interactions between oseltamivir and residues E119, D151, R152, R292, N347, and R371 across each simulation run. (B) through (D) are schematic views depicting the orientation of protein sidechains which form protein-drug hydrogen bonds.

Beyond disrupting the drug-protein hydrogen-bonding network, another mechanism through which protein mutations may induce drug resistance is by disruption of the hydrophobic packing of the drug into the protein binding pocket. Through inspection of the static crystal structures of the H274Y and N294S mutants of H5N1, it has been speculated that the mutations disrupt favourable hydrophobic packing interactions necessary for strong binding of oseltamivir. In our wild type simulations for both the H5H1 and H1N1pdm systems, the packing of oseltamivir’s pentyl moiety tended to favour close association with residues I222, R224, A246, and E276. To test the effect of mutations H274Y and N294S on hydrophobic interactions of oseltamivir’s pentyl group with the proteins, we calculated the solvent accessible surface area (SASA) of oseltamivir’s pentyl group for all simulation trajectories. While there was no significant change to the pentyl group SASA (henceforth referred to as PG-SASA) in the wild type and N294S mutant for either H5N1 or H1N1pdm neuraminidases, an outward rotation of the pentyl group was observed in the H274Y mutant simulations, visible in oseltamivir’s binding pose and evident in a notably higher calculated PG-SASA. The PG-SASAs for all simulated systems are plotted in Figure 5, with inset images of oseltamivir’s binding pose in simEQ2 (Figure 5A) and simEQ4 (Figure 5B) illustrating the rotation of the pentyl group towards the open mouth of the binding pocket.

Figure 5.

The solvent accessible surface area of oseltamivir's pentyl group (PG-SASA) in H5N1 and H1N1pdm wild type and mutant simulations. The PG-SASA in the N294S mutant simulation (simEQ5 and simEQ6) did not vary significantly from wild type simulations (simEQ1 and simEQ2). However, there was a huge loss of hydrophobic interaction between drug and protein in the case of the H274Y mutant simulations (simEQ3 and simEQ4), resulting in an increase in measured PG-SASA.

Previously published MD simulations performed over relatively short time scales (3 to 6 ns) have suggested two possible mechanisms: 1) that the H247Y mutation reduces the size of the hydrophobic pocket within the SA binding pocket near Oseltamivir’s pentyl moiety (Malaisree et al., 2008), and 2) that the H274Y mutation breaks a critical salt bridge between E276 and R224 to disrupt drug binding (Wang et al., 2009). Our longer (40ns) simulations were able to corroborate the former suggested mechanism (as shown by the increase in PG-SASA in the case of the H274Y mutant simulations) but not the latter. In fact, in all six of our simulations, E276 maintains stable charge-charge interactions (salt bridging) with R224 despite displacement of the drug from the protein I222-R224-A246-E276 pocket in the case of H274Y mutants. This drug displacement increases water penetration into the I222-R224-A246-E276 pocket (Figure 5B). Evidence from our simulations therefore supports predictions from earlier studies of a possible mechanism for the H274Y mutation-induced drug resistance through water infiltration and destabilization of favourable drug packing. However, because no change in PG-SASA was observed during our simulations in the case of N294S for either protein, N294S induced drug resistance is probably due to a different mechanism.

Figure 6.

Electrostatic surface potential of avian H5N1 and swine H1N1pdm neuraminidases in oseltamivir-bound simulations, revealing a positively charged pathway into the binding pocket. (A) through (C) show oseltamivir bound to H5N1 wild type, H274Y, and N294S drug-resistant mutant structures, respectively. (D) through (F) show oseltamivir bound to wild type, H274Y, and N294S drug-resistant mutant structures, respectively. The outer columns show a close-up view of the binding pocket, highlighted as a subset of the enter protein shown in the central column. The positions of the mutant residues are shown in green for residue 274 and 294 in the mutant systems

While alterations in hydrophobic sidegroup packing may explain in part the mechanism behind H274Y mutation-induced drug resistance, they fail to shed light on the role that the N294S mutation plays for oseltamivir binding inhibition. Initial investigations of the electrostatic surface potentials of drug-bound N1 neuraminidases have proposed that drug binding affinity is closely related to favourable charge-charge interactions with the electrostatic potential of the binding pocket wall, which exhibits a weak negative charge. In order to understand whether charge-charge interactions may play a role in mutation induced drug resistance, the electrostatic potentials were calculated and mapped onto the surfaces of the proteins. Even though prior studies have investigated the possible role of the electrostatic surface potential for drug binding in neuraminidases, the extensive electrostatic calculations required for fully understanding role of this phenomenon over a long simulation trajectory has not yet been done. The electrostatic surface potentials of the equilibrated systems were calculated and averaged across every trajectory frame in six equilibrium simulations. Extensive electrostatic calculations for the six systems show that the mouth of the binding pocket is positively charged, except for a narrow pathway of negative surface charge that seems to direct a possible binding pathway through which the drug may access the binding pocket. In Figure 6, the electrostatic surface potentials calculated from all six simulations is shown for the SA binding pocket with oseltamivir bound.

Mutation H274Y in the H1N1pdm virus shares the same sources of oseltamivir resistance as H5N1, including loss of hydrogen bonds with R152 and reduction of hydrophobic interaction of oseltamivir’s pentyl group. This leads to our suggestion to replace the pentyl group by a more hydrophilic one by adding a few hydroxyl groups on this bulky hydrophobic group. Our results do not support previous suggestions that the N294S mutation of H5N1 destabilizes the hydrogen bonding between oseltamivir and Y347, or that it disrupts the hydrophobic pocket, leading to drug resistance. Sources of drug resistance in mutant N294S remain unclear. Quantitative analysis of drug binding to the pocket is surrounded by a highly positive potential ring (colored blue). Simulations revealed the presence of a negatively charged pathway at the mouth of the binding pocket which may play a role for drug binding and mutation induced resistance, as the position of residue 294 maps directly onto the pathway while position of residue 274 is positioned on or adjacent( Le et al., 2009). WT and the mutants will be importantly complementary to obtained results in this chapter. Since the two mutations are nonactive-site, endpoint interactions alone cannot account for all drug resistance. Kinetics of drug binding is therefore also important for drug resistant mechanism. Especially, based on the decrease in association rate constants (K_on) of oseltamivir with H5N1 neuraminidase mutants (2.52 μM−1s−1 in the WT, 0.24 2.52 μM−1s−1 in H274Y and 1.1 2.52 μM−1s−1 in N294S mutants), we speculate that mutations also prevent drugs from entering the SA binding site. It is hoped that our observation of a possible binding pathway for oseltamivir will encourage further investigations that test the hypothesis posed here, by identifying whether the actual drug binding pathway follows this route.

4. Identification of drug binding pathway

4.1. Characteristics of the electrostatic surface potential of N1 neuraminidases shed light for further study on drug binding and drug resistance

Up to this point, all proposed mechanisms for oseltamivir resistance have focused on the effects of the mutations on the SA binding site but little has been known about the effects of those mutations on the kinetics of drug binding. Since both H274Y and N294S are nonactive site mutations, prior studies which focused on endpoint interactions between drug and proteins have been unable to provide a full understanding of how these mutations directly impact drug binding. In fact, given that the drug binding kinetics of H5N1 mutants are significantly diminished, it is possible that these mutations alter the binding process, and not necessarily just the specific endpoint interactions. It is well known that electrostatic surface potential of a protein can be an important driving force directing the diffusion of ligands into a protein’s active site. The resulting electrostatic maps shown in Figure 7A for H5N1 and Figure 7B for H1N1pdm reveal a highly negatively charged column of residues that forms a pathway 10 Å in length between the SA binding site and the edge of the binding cavity mouth. Electrostatic calculations also reveal that oseltamivir has a highly positive electrostatic surface potential, illustrated in Figure 7C. The question arises whether the negatively charged surface column plays a role in the binding and unbinding of oseltamivir, given a possible mutual attraction between oseltamivir and this column. To answer this question, we employed SMD simulations (described in the Methods section) to pull oseltamivir out of the SA binding site and probe possible unbinding pathways. Such simulations have

Figure 7.

Electrostatic surface potential of the SA binding pocket of H1N1pdm and oseltamivir. Shown in A) and B) are closeup views of the SA binding pocket with drug bound H1N1pdm and avian H5N1 neuraminidase, respectively. The region of the binding pocket where the drug bound possesses a highly negative potential (colored red), whereas the opening of the pocket is surrounded by a highly positive potential ring (colored blue). In C), a detailed surface electrostatic potential for oseltamivir. Shown are the “front” side facing the annulus of the binding pocket, and “back” side facing the interior of the binding pocket is shown.

4.2. Computational details

Six systems were modelled and simulated for oseltamivir bound H5N1, and H1N1pdm wild type and H274Y and N294S mutants. In simSMD1, steered molecular dynamics (SMD) simulations were used to remove oseltamivir from its stable binding site in H5N1 neuraminidase. In simFEQ1-10, equilibration simulations used a starting point generated from simSMD1 in which Oseltamivir has undergone an axial rotation such that it is partially displaced from its binding site. In total, 680 ns of simulations were carried out on system sizes of approximately 35,000 atoms.

The “Ensemble” column lists the variables held constant during the simulations; N, p, T, and V correspond to the number of atoms, pressure, temperature, and volume respectively. Under “Type”, EQ denotes equilibration, and SCV denotes constant velocity SMD simulation with the speed of 0.5 Å/ns. SimSMD1 is a steered MD simulation with the starting structure from equilibrated simEQ1 (“preflip” drug position). In simFEQ1 to simFEQ10, the starting structure reflected a “flipped” position of oseltamivir from simSMD1 after 7.5ns of simulation, whereas several of its main stabilizing hydrogen bonds to the protein had already been ruptured. 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 23 and ionized by NaCl (0.152M) to mimic physiological conditions.

All simulations were performed using NAMD 2.7 and the CHARMM31 force field (Mackerrel et al., 1998). 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 Nose´-Hoover Langevin piston method with a decay period of 100 fs and a damping time constant of 50 fs. 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 with a grid size of less than 1Å, along with the pencil decomposition protocol where applicable. SMD simulations fixed the center of mass of neuraminidase -carbons and applied a force to the center of mass of oseltamivir, along a vector connecting the two centers of masses. In simSMD1, a constant velocity protocol was employed, with stretching velocity of 0.5Å/ns. For the SMD spring constant we chose k_SMD = 3kBT/Å2 which corresponds to an RMSD value of (pkBT/k_SMD)1/2 < 0.6 Å.

4.3. Results

To simulate the binding of the drug would be the most natural approach, but the computations would be rather expensive. It is possible to the put drug in front of the receptor active site and the attraction forces between the drug and the receptor will drive it to the binding site. This process usually requires microsecond simulations and there is a high chance that the drug will fluctuate around the water box. Unbinding simulations however are feasible and may also reveal features that are characteristic for binding. In simSMD1, a pulling force was applied to rupture all the stabilizing hydrogen bonds between H5N1 and oseltamivir, in order to draw the drug away from the SA binding site. The results of simSMD1 show that the behaviour of oseltamivir under effect of a pulling force can be divided into three distinct stages: 1) from 0 to 8 ns, a buildup of forces during which hydrogen bonds between oseltamivir with E119, D151 and R152 are destabilized; 2) at 8 ns where the remaining stable hydrogen bonds with R292 and R371 rupture; 3) after 8 ns during which the drug is pulled out of the binding pocket. Figure 8 shows the force dependent rupture of hydrogen bonds in simSMD1 by plotting both hydrogen bond stability and (in inset) the force vs. time curve. To our surprise, it was observed that

Figure 8.

Distances between hydrogen bond acceptor-donor pairs between oseltamivir and active site amino acids vs. simulation time in simSMD1. Most hydrogen bonds were quickly broken by the pulling force except for those with R292 and R371, which fully ruptured only after 8 ns, corresponding to the peak of the curve of the applied force vs. simulation time (shown in inset). Position of oseltamivir and all its possible hydrogen bonding pairs with residues located along this pathway, it was seen that nonspecific electrostatic attractions formed the predominant interactions between drug and protein.

Oseltamivir followed a lateral escape path through strong interaction with the negatively charged column of residues indentified via electrostatic mapping above, despite application of force to pull the drug straight out of the binding pocket. This divergent path taken by the drug is shown in Figure 9A-D. A second notable observation in simSMD1 was that just before rupture of hydrogen bonds between oselvamivir’s carboxyl functional group and residues R292 and R371 of the binding pocket, the drug underwent a rather significant rotation due to the destabilization of hydrogen bonds with E119, D151 and R152. It is likely that this rotation or “flip” is crucial for placing oseltamivir in a position which permits it to leave the SA binding site. Despite application of force to pull oseltamivir straight out of its binding pocket, the drug followed a somewhat lateral path through the electrostatically charged funnel. The “preflip” or stably bound state for oseltamivir in H5N1 is shown in Figure 10A1, 10A2, 10A3, and the flipped transition state shown in Figure 10B1 and 10B2, 10B3.

It was observed in simSMD1 that following the transition to the flipped state, very little force was then required to subsequently draw oseltamivir out of its binding pocket in neuramidase. This result suggests that one may be able to probe the unbinding pathway via diffusion, if the starting state consists of oseltamivir already in its flipped orientation. We therefore performed additional equilibrium simulations (simFEQ1-10) with oseltamivir already in this flipped state. From these simulations, we were able to observe two distinct outcomes: 1) oseltamivir is able to escape the binding pocket by interacting favourably with the charged binding funnel and 2) the drug returns, not unexpectedly, to its stably bound “preflip” state. Each simulation for FEQ1-10 was run long enough to observe either outcome, with the exception of FEQ5. SimFEQ5 was a special case in which the drug, after following the binding funnel to escape the protein, actually rebound to the SA active site through the same binding funnel after we extended the simulation to follow its movement. A summary of observed outcomes from these simulations is shown in Table 2. In simFEQ1-5, oseltamivir successfully escaped the SA binding site, whereas in simFEQ6-10, oseltamivir returned to its stably bound “preflip” state. Oseltamivir was observed to diffuse out of the SA active site after strong interaction with the electrostatically charged binding funnel (described above) in five out of ten equilibrium simulations we performed (simFEQ1-5).

Figure 9.

Forced unbinding of oseltamivir from H5N1 neuraminidase. Shown here are the relative positions of oseltamivir on the electrostatic surface of avian H5N1 neuraminidase during simSMD1.At 0 ns (A), oseltamivir was stably bound within the active site, as seen in simEQ1. Application of force ruptured the stabilizing hydrogen bonds between H5N1 and oseltamivir (see Figure 9), drawing the drug away from its stable binding site within 10 ns, as shown in B. Over the next 2.5 ns of pulling, oseltamivir followed the charged binding funnel (shown in C) until it was completely free of the protein binding pocket after 15 ns, as shown in D).

In four of the cases (simFEQ1-3, 5) oseltamivir was observed to diffuse along the full length of the binding funnel before disassociating with neuraminidase. In simFEQ5 we observed not only a diffusion of oseltamivir through the charged binding funnel, but also the re-entry.

Snapshots from each of these events in FEQ5 are shown in Figure 10. Analysis of interactions of the newly rebound oseltamivir with active site residues from 50 to 100 ns revealed that the drug was stabilized by hydrogen bonds with Y406, R292, D151, E119, and R118, even though the pentyl group had not yet moved to its requisite hydrophobic pocket (I222-R224- A246-E276).

Figure 10.

Oseltamivir in “flipped” position (position of oseltamivir at 7.5ns in simSMD1) in comparison with its stable equilibrium position,

Name	Result	Time (ns)
simFEQ1	Drug escape via binding funnel	15
simFEQ2	Drug escape via binding funnel	10
simFEQ3	Drug escape via binding funnel	15
simFEQ4	Drug interaction with binding funnel but escape via 430-cavity	50
simFEQ5	Drug escape and rebind into SA pocket via binding funnel	100
simFEQ6	Drug returned to “preflip” position	10
simFEQ7	Drug returned to “preflip” position	15
simFEQ8	Drug returned to “preflip” position	50
simFEQ9	Drug returned to “preflip” position	50
simFEQ10	Drug returned to “preflip” position	50

Table 2.

Summary of FEQ1-10 simulations starting from “flipped” position of oseltamivir taken from simSMD1 at 7.5 ns.

Figure 11.

Escape and rebinding of oseltamivir through the electrostatic binding funnel in H5N1 neuraminidase during simFEQ5. Shown here are snapshots of simFEQ5, in which oseltamivir first diffused out of the SA active site through interaction with the electrostatic binding funnel within the first 25ns of simulation. Between 28 and 35ns, oseltamivir diffuses and approaches the periphery of the binding pocket away from the binding funnel, but is prohibited from entering due to electrostatic repulsion (45ns). However, between 45 and 50ns, oseltamivir was observed to approach and enter the neuraminidase binding pocket through the binding funnel, adopting a stable position within the SA binding pocket through hydrogen bonds with Y406, R292, D151, E119, and R118.

The result has shed light on the important role of the electrostatic surface potentials in directing the diffusion of oseltamivir into the SA binding site. From our simulations, it is clear that the negatively charged funnel serves as a prominent binding and unbinding pathway for oseltamivir in the wild type systems investigated with SMD and followup equilibrium simulations in simSMD1 and simFEQ1-10. It turns out also that the binding funnel may possibly play a crucial role in drug resistance caused by mutations. The conspicuous location of residue 294, which maps directly onto this negatively charged pathway, may play a key role in the N294S mutation for disrupting the proper guidance of the drug into its binding pocket. Thus, it is possible that the 274 and 294 mutations may confer drug resistance by not only disrupting the end-point interactions of oseltamivir but also its entry into the binding pocket by interfering with the binding funnel. While sources of end-point interactions, including hydrogen bonds, hydrophobic packing and solvent infiltration, of oseltamivir resistance have been thoroughly studied, little is known about the kinetics of the drug binding in mutants at atomic level. The idea that drug resistant mutants actually disrupt entry of oseltamivir into the SA active site of neuraminidase through disruption of an electrostatic binding funnel is in part supported by experiments which have noted reduced drug binding kinetics in H5N1 H274Y and N294S mutants. Even though the oseltamivir-resistant mutations were seen located in or adjacent to the funnel, clearly additional study is still needed for a full understanding of how the H274Y and N294S mutations weaken the binding of the drug. Furthermore, the mutations might not only affect the electrostatic gradient but also the geometry of the drug binding and unbinding path. Future or followup studies should therefore focus on the specific drug entry/exit pathways for oseltamivir-resistant mutant systems, sampling timescales great enough (such as in the case of simFEQ5 described above) to capture both the binding.

5. Challenges and opportunities

We have presented three important applications of MD simulation in drug design against influenza A virus. While MD simulation has been proven as a powerful tool complement to experiment in drug discovery processes, the problems in accuracy of drug parameter and limitation in time scale of MD simulation definitely affect applications of this tool. The challenge in drug force field development is that, unlike uniformed amino acids or nucleotides, each drug has its own structure and requires separate parameter. Since force field for general macromolecules have been well developed, attention should be shifted to development of adequately accurate drug force field (Wang et al., 2004; Zoete et al., 2012). Another challenge in MD simulation is limitation of time scale. Several important phenomena which are important for drug design application happen in the time scales that can’t be done by even the biggest supercomputers in the world. Over the past decade several methods were developed to speed up classical MD simulation including SMD, accelerated MD (Hamelberg et al., 2009), coarse-grained MD simulation (Izvekov et al., 2005). Last but not least, the implementations of specialized hardware including CBE, GPU, and FPGA architectures are important to broaden application of MD simulation to modern drug discovery (Horacio et al., 2012). The recent finding of oseltamivir binding pathway using GPU acceleration and SMD (Le et al., 2010) is an example that the combined improvement of theoretical algorithm and specialized hardware architecture has enabled us to increase the success rate and decrease the cost in rational drug discovery.

Acknowledgement

The author wish to thank the Institute for Computational Science and Technology at Ho Chi Minh City and NAFOSTED (The National Foundation for Science and Technology Development) for partially support the research presented in this chapter.