Selective Binding of Antiinfluenza Drugs and Their Analogues to

Sep 22, 2010 - intrinsically lower in energy than the closed conformation and that oseltamivir (tamiflu) favors binding to the closed conformation thr...
3 downloads 0 Views 5MB Size
12958

J. Phys. Chem. B 2010, 114, 12958–12964

Selective Binding of Antiinfluenza Drugs and Their Analogues to ‘Open’ and ‘Closed’ Conformations of H5N1 Neuraminidase Pei Wang§ and John Z. H. Zhang*,†,‡,§ Institute of Theoretical and Computational Science, and State Key Laboratory of Precision Spectroscopy, Department of Physics, East China Normal UniVersity, Shanghai 200062, China, and Department of Chemistry, New York UniVersity, New York, New York 10003 ReceiVed: April 4, 2010; ReVised Manuscript ReceiVed: September 3, 2010

It was suggested that the open conformation of the 150-loop of H5N1 avian influenza neuraminidase is intrinsically lower in energy than the closed conformation and that oseltamivir (tamiflu) favors binding to the closed conformation through a relatively slow conformational change [Russell, R. J. Nature 2006, 443, 45-49]. In the present work, a systematic computational study is performed to investigate the binding mechanism of five ligands to H5N1 neuraminidase (H5N1 NA) with the 150-loop in both open and closed conformations through molecular docking, molecular dynamics simulations, and MM/PBSA free energy calculation. Our result shows that the electrostatic interactions between polar groups on the 150-loop and the charged groups of the ligands play a key role on the binding selectivity. In particular, ligands having a small positively charged group favor binding to the closed conformation of H5N1 NA, while those having a large positively charged group generally prefer binding to the open conformation. Our analysis suggests that it may be possible to design new inhibitors with large basic groups that are selective for the open conformation and thereby have stronger binding affinity to H5N1 neuraminidase. 1. Introduction Avian influenza virus can cause influenza in birds and some mammals. Subtypes of influenza virus are labeled according to two functional membrane proteins: hemagglutinin (HA) and neuraminidase (NA). The former allows the recognition of target cells through binding to cells’ sialic acid-containing receptors. The latter cleaves terminal sialic acid residues from carbohydrate moieties on the surfaces of infected cells and thereby promotes the efficiency of virus release from cells.1 H5N1, one of the most dangerous subtypes of the influenza A virus which can cause illness in humans and many other animal species, has became one of the world’s major influenza pandemic threats according to the reports of the World Health Organization (WHO).2 Therefore, the study of inhibitors or vaccines of H5N1 has attracted interest and worldwide attention. Neuraminidase (NA) is a popular drug design target. Through structure-based enzyme inhibitor design, two antiviral drugs, oseltamivir (Tamiflu)3 and zanamivir (Relenza),4 have been produced and used as clinical treatments for influenza. In addition, peramivir, an experimental antiviral drug, has also advanced into phase II studies.5 All of these inhibitors act as a transition-state analogue inhibitor of influenza neuraminidase and thereby prevent new viruses from emerging from infected cells. These successful drugs have helped to discover the active sites of various kinds of neuraminidases and reveal that these enzymes are comparatively rigid with only insignificant conformational changes of the sites upon inhibitor binding. Consequently, the active sites of all influenza neuraminidases * To whom correspondence should be addressed. E-mail: john.zhang@ nyu.edu. † Institute of Theoretical and Computational Science, East China Normal University. ‡ State Key Laboratory of Precision Spectroscopy, Department of Physics, East China Normal University. § Department of Chemistry, New York University.

contain three arginine residues, Arg 371, Arg 292, and Arg 118, that orient and bind to the carboxyl group of sialic acid, as well as Arg 152 that forms a hydrogen bond with acetylamino group of the substrate.6 Russell et al. conjectured that a 150-loop (residues 147-152) (Figure 1) of apo-H5N1 NA is more stable in the ‘open’ conformation but that the ‘closed” conformation is preferred when oseltamivir is bound to H5N1 NA.7 However, it is not clear how oseltamivir binds to these two conformations and what the respective binding affinity is. Thus, taking into account the binding selectivity for the open or closed conformation of the 150-loop is important for new antiinfluenza drug design. Currently, experimental methods cannot provide information on the difference in binding affinity of ligands to the open and closed conformations of H5N1 NA. Also, various recognition elements of different inhibitors are not apparent, and their inhibitory activities are roughly equal.8 To explore the binding properties of H5N1 NA for the purpose of inhibitor design, it is desirable to perform a detailed analysis of how open and closed conformations of H5N1 NA recognize different ligands. In this study, we carry out systematic computational studies to investigate possible ligand-H5N1 NA binding complexes through molecular docking, molecular dynamics (MD) simulations, and molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) free energy calculations.9,10 In addition to three NA inhibitors mentioned above, we also include two additional analogues of oseltamivir with inhibitory activity, tamiphosphor and 13b, which was synthesized by Wong’s group.11 The chemical structures of ligands in this study are shown in Figure 2. The purpose of our study is to determine the affinity differential of different ligands when binding to the open and closed conformations of H5N1 NA. The information derived from such analysis could enable discovery of new antiinfluenza drugs, especially those with better binding selectivity to specific conformations of H5N1 NA.

10.1021/jp1030224  2010 American Chemical Society Published on Web 09/22/2010

Selective Binding of Antiinfluenza Drugs

J. Phys. Chem. B, Vol. 114, No. 40, 2010 12959

Figure 1. (A) A superposition of open (2HU0: blue) and closed (2HU4: yellow) conformations of H5N1 NA. (B) A view of oseltamivir (red) in the active cavity of 2HU0 with the adjacent 150-loop. (C) A view of oseltamivir (red) in the active cavity of 2HU4.

Figure 2. Chemical structures of ligands in this study.

2. Computational Details Atomic coordinates of H5N1 neuraminidase were obtained from the Protein Data Bank (PDB): open (2HU0) and closed (2HU4) conformations of H5N1 NA. Protonation states for histidine residues were determined at an apparent pH 6.5 using the PDB2PQR web server.12 Each of the five ligands was first optimized with Gaussian03 at HF/6-31G* level of theory. Single-point calculations were then performed to obtain the

electrostatic potential around each compound using the same basis set and level of theory as in the optimizations step. The atomic charges of the ligands used for molecular mechanics calculations were derived from the ESP by using the RESP13 program implemented in the AMBER 9.0 package. In view of the crystallographic conditions of H5N1 NA, all ionizable groups both on the ligand and protein were configured in their characteristic ionized states at pH 6.5. Oseltamivir, zanamivir,

12960

J. Phys. Chem. B, Vol. 114, No. 40, 2010

and peramivir were given a formal charge of 0 while tamiphosphor and 13b had a charge of -1. The protein chains bound to oseltamivir in two original crystal structures (chain B in structures 2HU0 and 2HU4) were used as receptors for different small ligands. AutodockTools version 1.5.0 was used to add polar hydrogen and assign Gasteiger charges14 to each of these two NA structures. AutoGrid version 4.0 was used to create affinity grids centered on the active site. Each grid enclosed an area of 78 Å by 66 Å by 56 Å with 0.375 Å spacing, and affinity grids were calculated for all the following atom types: A (aromatic carbon), C, N, NA (hydrogenbond-accepting N), O, OA (hydrogen-bond-accepting O), P, S, SA (hydrogen-bond-accepting S), H and e (electrostatic). Autodock version 4.0.1 with the Lamarckian genetic algorithm was used to simulate ligand-receptor docking.15 Different binding modes with the best scores sorted by the lowest energy and the largest populations for each molecule were selected for further MD study and analysis. Molecular dynamics simulations were carried out using the SANDER module of the AMBER9.0 program with the ff99SB version of Amber force field.16 Generalized AMBER force field (GAFF) provided parameters for small organic molecules to facilitate simulations of ligands in conjunction with biomolecules, and we also replaced those parameters of a guanidinium moiety with those generated from arginine. For each selected binding mode obtained from AutoDock, about 35 000 TIP3P water molecules with a 10 Å buffer were added around the complex. Four or five Na+ counterions were added to maintain the neutrality of the system. Production MD was performed for 3 ns with a 2 fs time step, and the SHAKE algorithm17 was employed to fix bond lengths involving hydrogen atoms. The particle-mesh Ewald method (PME)18 was used to provide an efficient and effective treatment of long-range electrostatics interactions. During the simulation, nonbonded cutoffs were set to 10 Å; the temperature was maintained at 300 K. After minimization of 1500 steps and equilibration for 50 ps with the system gradually warmed to 300 K, complex conformations were collected every 1 ps for subsequent 3 ns simulations. The binding free energies of all the systems as well as component energy contributions were calculated using a molecular mechanics/Poisson-Boltzmann surface area (MM/ PBSA) approach.10 The average free energy of the complex, receptor, or ligand is composed of the molecular mechanic, solvation, and entropic energies of the system over the trajectory. In this study, the single trajectory approach applied19 as an estimation of energies in this manner has proven successful in many studies.20,21 This approach involves the extraction of thermodynamic data from a single trajectory of the protein-ligand complex based on the assumption that there is no significant conformational change between bound and unbound state and that error calculation in this manner can somewhat hide the effect of incomplete sampling. The snapshots of each system were sampled from the last 1 ns single trajectory with an interval of 10 ps, i.e., 100 snapshots. Molecular mechanics energy was obtained using the Amber 9 Sander module, and the electrostatic contribution to the solvation energy was calculated with the PBSA program of AMBER 9.0. Nonpolar contributions to solvation were estimated as a function of the solvent-accessible surface area (SASA). From those simulations, the nonpolar part of the solvation energy was estimated using the simple empirical relation, where the solventaccessible surface area (SASA) is estimated using Sanner’s algorithm implemented in the MSMS program22 with a solvent probe radius of 1.4 Å and the PARSE atomic radii parameters.23

Wang and Zhang The values for γ and b were obtained experimentally from the transfer of small hydrocarbons24 and are 0.00542 kcal/mol · Å2 and 0.92 kcal/mol, respectively. One normal mode calculation was performed by using Nmode module in the AMBER 9.0 package. For each complex, six snapshots from the MD trajectories were taken at 200 ps intervals. The normal mode calculation to estimate the entropy contribution is somewhat problematic and time-consuming25,26 as is the reference part of binding free energy estimation. 3. Results and Discussion Docking of oseltamivir and its four analogues to the active site of both open and closed conformations, respectively, of H5N1 NA was performed. In general, molecular docking gives quite different lowest energy orientations, especially for a molecule with many degrees-of-freedom (DOF) and weak binding. Thus, cluster analysis was performed using a tolerance of 2.0 Å rmsd (root-mean-square deviation) in order to identify the most favorable docking modes from the docking results. However, in this study, almost all the docking modes were statistically significant (see Supporting Information Table S1). The optimal binding modes of each molecule with the best scores and largest populations were selected for further MD study, and analysis focused on the active site, the SA cavity,27 adjacent to the 150-cavity.7 Based on the initial docking result, a detailed MD analysis of five ligand-H5N1 NA complexes was performed and their binding free energies were calculated for the binding, respectively, to the open and closed conformations of H5N1 NA. To obtain more accurate information on the binding mode and ligand-protein interactions, a series of 3 ns MD simulations were carried out. Figure 3 plots the rmsd of the protein backbone atoms of 2HU0 and 2HU4 in complex with five ligands. For each complex, the optimal docking modes of ligands were used as the initial structure. Interestingly, all the complexes have similar rmsd values and reached equilibrium after 0.5 ns simulation. Oseltamivir and Tamiphosphor Bind to the Closed Conformation of H5N1 NA. To obtain information about the selectivity of binding of these five ligands to the open and closed conformations of H5N1 NA, the MM/PBSA method was used to calculate the binding free energies. Table 1 summarizes the estimated free binding energies along with the contribution of individual energy terms for each ligand-protein complex. Our free energy calculation shows that oseltamivir prefers binding to the closed conformation with a binding free energy of -8.5 kcal/mol. For comparison, the binding free energy of oseltamivir to the open conformation of H5N1 NA is +12.4 kcal/mol. This significant difference in binding free energy is almost entirely due to the difference in electrostatic interaction of oseltamivir in different conformations. Our calculation thus confirmed the suggestion based on the experimental crystal structure data that oseltamivir favors binding to H5N1 NA in closed conformation, which is more energetically favorable. This significant electrostatic change also suggests that the conformational change of the 150-loop upon binding of oseltamivir is mainly driven by electrostatic interactions at the binding site, which will be discussed later. Our calculation also finds that tamiphosphor, the phosphate analogue of oseltamivir, shows some preference of binding to the closed conformation of the 150-loop with a difference in binding free energy of ∼4 kcal/ mol. Analysis of binding interaction indicates that tamiphosphor binds to the closed conformation with a more favorable van der Waals (vdW) interaction by ∼8 kcal/mol over that of the

Selective Binding of Antiinfluenza Drugs

J. Phys. Chem. B, Vol. 114, No. 40, 2010 12961

Figure 3. The rmsd of the protein backbone atoms (C, CR, and N atoms) of the complexes of (A) H5N1 NA in open and (B) closed conformation and the set of five ligands: oseltamivir, zanamivir, peramivir, tamiphosphor, and 13b.

TABLE 1: Various Energy Contributions to the Free Energy of Binding for Five Ligands to H5N1 NAa ligand oseltamivir zanamivir peramivir tamiphosphor 13b

protein

∆Eelecb

∆EVDWc

elec d ∆Gsolv

nonpolar e ∆Gsolv

elec f ∆Eint+solv

∆E

T∆Sg

∆Gh

2HU0 2HU4 2HU0 2HU4 2HU0 2HU4 2HU0 2HU4 2HU0 2HU4

-141.9 -167.6 -199.1 -214.7 -197.5 -152.8 -205.6 -213.4 -260.8 -269.1

-25.4 -27.1 -27.4 -24.6 -27.6 -27.3 -20.1 -28.1 -22.5 -24.5

157.4 163.6 191.2 203.5 175.4 140.5 184.8 201.2 229.9 276.1

-4.8 -4.6 -4.6 -4.6 -4.9 -4.8 -4.6 -4.9 -4.8 -4.6

20.2 0.8 -3.1 -6.5 -17.3 -7.6 -16.2 -7.3 -26.1 11.6

-9.9 -30.9 -35.2 -35.7 -49.8 -40.0 -41.0 -40.3 -53.4 -17.5

-22.3 -22.4 -21.9 -22.7 -21.7 -21.5 -28.9 -24.2 -24.3 -23.8

12.4 -8.5 -13.3 -13.0 -28.1 -18.5 -12.1 -16.1 -29.1 6.3

a

All energies are in kcal/mol. b Nonbonded electrostatics. c Nonbonded van der Waals. d Electrostatic component to solvation. e Nonpolar component to solvation. f Total electrostatic change upon binding. g Entropic contributions to binding. h Total change of free energy in binding.

open conformation. To state it in another way, the relatively small size of tamiphosphor enables it to fit better to the closed conformation of H5N1 NA. Zanamivir, Peramivir, and 13b Prefer the Open Conformation of H5N1 NA. We also find that zanamivir, peramivir, and 13b have significant selectivity for the open conformation of H5N1 NA. In detail, the overall electrostatic and vdW components to binding provide thermodynamic benefit to these ligands. The nonbonded electrostatic potential of zanamivir binding to the closed conformation of H5N1 NA is 15.6 kcal/ mol lower that that of the binding to the open conformation of H5N1 NA. It appears that zanamivir gains this electrostatic preference by having a large positively charged guanidinium group that interacts favorably with Asp 151 on the 150-loop and Glu 229 at the bottom of the active site. But the electrostatic component to solvation leads to a loss of 12.3 kcal/mol for zanamivir binding to the closed conformation of H5N1 NA. On the other hand, zanamivir binding to the open conformation of H5N1 NA has a more beneficial vdW component by ∼3 kcal/mol over that of the closed conformation. Thus zanamivir seems to have considerable binding selectivity to the open conformation of H5N1 NA. Our results also reveal that peramivir and 13b selectively bind to the open conformation of H5N1 NA. The electrostatic energy of peramivir decreases by ∼45 kcal/mol upon binding to the open conformation from the closed conformation. Although this decrease is partially offset by the decrease of ∼35 kcal/mol from the electrostatic component of solvation, the free energy affinity difference for peramivir between binding to open and closed H5N1 NA conformations is up to 9.6 kcal/mol. On the contrary, the nonbonded electrostatic energy of 13b increases by ∼9 kcal/ mol upon binding to the closed conformation from the open

conformation because of strong interactions between carboxylate-guanidinium and phophonate-guanidinium ion pairs, but this increase is overwhelmed by the growing electrostatic component of solvation at around 46 kcal/mol. These differences between nonbonded electrostatic energy and solvation electrostatic energy are most likely due to the lack of flexibility of zanamivir, peramivir, and 13b that precludes optimal charge complementarities in the active site to offset desolvation penalties when binding to the closed conformation of H5N1 NA. More interestingly, these three ligands, which prefer the open conformation of H5N1 NA, have the large, positively charged guanidinium group. Furthermore, according to the vdW interaction change analysis, it is believed that a large group is more favorable to fit into the open conformation of H5N1 NA, which has a relatively large binding site. Antiinfluenza Effect of Phosphonate Group of Tamiphosphor and 13b. The phosphate group is generally applied as a bioisotere of carboxylate in drug design.28 Our results agree with the experimental finding that these two phosphonate analogues, tamiphosphor and 13b, of oseltamivir are significantly more potent than the carboxylate congeners against H5N1 NA. Indeed, the estimated ∆G obtained by using the relation ∆G ) -RT ln IC50 could be used for making a rough estimation of the free energy, because the IC50 has been shown to correlate with the constant of inhibition (Ki) for small inhibitors.29 All ligands are correctly rank-ordered and in comparable agreement with experimental data as shown in Table 2. Using this combined docking and MM/PBSA approach, we also predicted the binding selectivity of tamiphosphor and 13b; the former favors binding to the closed conformation of H5N1 NA while the latter tends to bind to the open one.

12962

J. Phys. Chem. B, Vol. 114, No. 40, 2010

Wang and Zhang

TABLE 2: Favorable Binding Free Energies with Experimental Valuesa

A view of the active site binding of 13b to H5N1 NA in the open 150-loop conformation is shown in Figure 4B. Like that of tamiphosphor, 13b is stabilized by the salt-bridge network with phophonate-Arg 371-Arg 292. The average distances of the phophonate-Arg 371 and phophonate-Arg 292 ion pairs are 2.76 Å ( 0.15 and 2.71 Å ( 0.08, respectively. Interestingly, because of the large, positively charged guanidinium group, 13b can interact with the 150-loop even in open conformation. Arg 152 of 2HU0 forms a hydrogen bond with the oxygen of the N-methylacetamide group of 13b with occupancy up to 99% along a 3 ns MD simulation. But this hydrogen bond is eliminated when 13b binds to 2HU4, the closed conformation of H5N1. Therefore, we can speculate that interaction with the flexible 150-loop mainly determines the binding selectivity for different ligands. It seems that ligands having larger, negatively charged groups can bind more strongly to the open conformation of H5N1 NA due to strong electrostatic interaction with the 150-loop, which leads to more optimal charge complementarities in the active site. Combining the free energy data and analysis of binding mechanism, we conclude that electrostatic contribution is a dominant factor in determining binding selectivity. Optimizing the oseltamivir structure with a more electronegative polar group may result in stronger binding. However, the high polarity of the phophonate and guanidinium group may interfere with oral bioavailability.11 Further investigation of formulations and the use of prodrugs is needed for further drug development. Hydrogen Bonding between Antiinfluenza Drugs and the 150-Loop of H5N1 NA. The ‘150-cavity’ is adjacent to the active site when the 150-loop is in open conformation7(Figure 1) but will shrink if the 150-loop is in closed conformation. Thus, the flexibility of the 150-loop of H5N1 NA (residues 147-152; Gly-Thr-Val-Lys-Asp-Arg) could influence the topology of the binding site and further affect the binding selectivity of ligand interactions. Upon superposition of the open and closed 150-loop conformations, as shown in Figure 1, the 150-loop tends to move to the active site in the closed conformation. Distance between Val 149 in these two conformations is about 7 Å, and the carboxyl group of Asp 151 moves toward to the active site. In order to investigate the influence of the 150-loop to binding selectivity, we examined electrostatic contributions of the 150-loop to binding for three antiinfluenza drugs when they interact with H5N1 NA. Figure 5 shows percent occupancies of the hydrogen bonds between the guanidinium group (nitrogen as the acceptor) of Arg 152 and the N-methylacetamide group (oxygen as the donor) of the ligands. Thus, the binding of osetamivir favors the conformational change of the 150-loop, which can help it form a hydrogen bond with Arg 152 and further improve its binding affinity. However, this movement of the 150-loop is disfavored by the binding of zanamivir and peramivir. It seems that because of the larger size of the guanidinium group compared to the amino group, these two ligands can interact with the 150-loop even in open conformation, while the hydrogen bonds are energetically unfavorable and are precluded upon binding to H5N1 NA with the 150-loop in a closed conformation. In addition, the conformational change of the 150-loop allows a salt bridge to form between Asp 151 and oseltamivir. The average distance between the carboxyl group of Asp 151 and the amino group of oseltamivir is 6.41 Å ( 0.76 when binding to the open 150-loop, whereas the corresponding average distance is 4.64 Å ( 0.35 when binding to the closed 150-loop. But the same

∆Gestimated ∆Gexperimentd

oseltamivirb

tamiphosphorb

13bc

-8.5 -9.8

-16.1 -10.8

-29.1 -12.0

a All energies are in kcal/mol. b Estimated binding free energy at the closed 150-loop conformation. c Estimated binding free energy at the open 150-loop conformation. d Experimental data given as IC50 for all ligands. For direct comparison to calculated affinities, conversion is made according to ∆G ) -RT ln IC50.

Figure 4. (A) A view of the binding mode of tamiphosphor in 2HU4 after a 3 ns MD simulation. (B) A view of the binding mode of 13b in 2HU0 after a 3 ns MD simulation. The dotted lines represent salt bridges and hydrogen bonds.

Figure 4A depicts a view of the active site binding of tamiphosphor to H5N1 NA in a closed conformation after a 3 ns MD simulation. Tamiphosphor was integrated into the salt bridge networks with phophonate-Arg 371-Arg 292 and aminoAsp 151-Glu 119-Glu 277, the negatively charged phophonate group and positively charged amino group playing central roles by forming salt bridges to stabilize the binding interactions. The average distances between phophonate-Arg 371 and phophonateArg 292 are 2.74 Å ( 0.08 and 2.76 Å ( 0.10, respectively. The average distances between amino-Asp 151, amino-Glu 119, and amino-Glu 277 are 2.92 Å ( 0.15, 3.01 Å ( 0.34, and 4.78 Å ( 0.23, respectively. Glu 277 also forms hydrogen bonds with the nitrogen atom of the amino group of tamiphosphor, with occupancy up to 95% along a 3 ns MD simulation. However, this occupancy value is just 75% during an MD simulation of tamiphosphor binding to the open conformation of H5N1 NA. It seems that Glu 277 is of a key importance for the binding selectivity of tamiphosphor.

Selective Binding of Antiinfluenza Drugs

J. Phys. Chem. B, Vol. 114, No. 40, 2010 12963 Acknowledgment. John Zhang acknowledges partial financial support from NYU Research Challenge Fund. Supporting Information Available: Table S1. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes

Figure 5. The percent occupancies of the hydrogen bonds between the guanidinium group (nitrogen as the acceptor) of Arg 152 and the N-methylacetamide group (oxygen as the donor) of the ligands. All hydrogen bonds are measured with 3.0 Å and 120° cutoffs for the distance and angle, respectively, and the mean values are obtained from 100 snapshots evenly extracted from the last 1 ns MD simulations, the same snapshots as those used for the MM/PBSA calculations.

improvement of binding is not observed for zanamivir and peramivir. Moreover, Asp 151 could form hydrogen bonds with the guanidinium group of zanamivir when binding to the open conformation of H5N1 NA, with occupancy up to 98% during simulation, while this occupancy value decreased to 12% when binding to the closed conformation. In general, the observation of the open conformation of H5N1 NA suggests that it is intrinsically lower in energy than the closed conformation for this enzyme. Thus, our work suggests that it is possible to design new inhibitors for H5N1 NA having a large, positively charged group, which are selective for the open 150-loop conformation and would thereby have the potential to bind more strongly than oseltamivir. Conclusions The binding complexes of oseltamivir and its four analogues with open and closed conformations of H5N1 NA were simulated from docking and MD calculations. Free energy analysis of the protein-ligand ensemble using the MM/PBSA approach was used to identify the most possible binding modes in order to study the selectivity differential. Our simulation shows that the conformation of the 150-loop is the key factor that influences the binding potential of the inhibitors. The open or closed conformational state can either change the size of the active site or affect the electrostatic interactions with ligands. In particular, in the open conformation, Asp 151 can form a hydrogen bond with a ligand possessing a large basic group, while in the closed conformation, this interaction is more favorable for oseltamivir binding. Combining the free energy calculation data and binding conformation analysis, we conclude that the phosphonate group can strongly contribute electrostatic interactions of binding due to extra negative charge. Overall, ligands with a small basic group, such as amino, favor the closed conformation of H5N1 NA. But for those inhibitors possessing a large, positively charged group, such as guandinium, binding to the open conformation of H5N1 NA is favored. As H5N1 undergoes mutation and reassortment, variations are created that can enable infection of species not previously known to carry the virus. Development of new antiflu drugs targeted to H5N1 and its mutations is necessary. While this study suggests binding modes of H5N1 NA in different conformations, this approach can be used as an analytical tool to design new drugs and predict the activity without crystal structure determination.

(1) Moscona, A. Drug Therapy: Neuraminidase Inhibitors for Influenza. N. Engl. J. Med. 2005, 353, 1363–1373. (2) World Health Organization. H5N1 avian influenza: timeline of major events, Dec 2008. (3) Kim, C. U.; et al. Influenza neuraminidase inhibitors possessing a novel hydrophobic interaction in the enzyme active site: design, synthesis, and structural analysis of carbocyclic sialic acid analogues with potent antiinfluenza activity. J. Am. Chem. Soc. 1997, 119, 681– 690. (4) von Itzstein, M.; et al. Rational design of potent sialidase-based inhibitors of influenza virus replication. Nature 1993, 363, 418–423. (5) National Institutes of Health. Evaluation of the efficacy and safety of peramivir in adults with acute serious or potentially life-threatening influenza, Dec 2008. (6) Clercq, E. Antiviral agents active against influenza A viruses. Nat. ReV. Drug DiscoVery 2006, 5, 1015–1025. (7) Russell, R. J.; et al. The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design. Nature 2006, 443, 45–49. (8) Hurt, A. C.; et al. Susceptibility of highly pathogenic A(H5N1) avian influenza viruses to the neuraminidase inhibitors and adamantanes. AntiViral Res. 2007, 73, 228–231. (9) Kuhn, B.; et al. Validation and use of the MM-PBSA approach for drug discovery. J. Med. Chem. 2005, 48, 4040–4048. (10) Masukawa, K. M.; Kollman, P. A.; Kuntz, I. D. Investigation of neuraminidase-substrate recognition using molecular dynamics and free energy calculations. J. Med. Chem. 2003, 46, 5628–5637. (11) Shie, J.; et al. Synthesis of tamiflu and its phosphonate congeners possessing potent anti-influenza activity. J. Am. Chem. Soc. 2007, 129, 11892–11893. (12) Dolinsky, T. J.; et al. PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res. 2007, 35, W522-5. (13) Bayley, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A wellbehaved electrostatic potential based method using charge restraints for deriving atomic charges: The RESP method. J. Phys. Chem. 1993, 97, 10269–10280. (14) Gasteiger, J.; Marsili, M. Iterative partial equalization of orbital electronegativitysa rapid access to atomic charges. Tetrahedron 1980, 36, 3219–3228. (15) Morris, G. M.; et al. Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. J. Comput. Chem. 1998, 19, 1639–1662. (16) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of multiple amber force fields and development of improved protein backbone parameters. Proteins: Struct., Funct., Genet. 2006, 65 (3), 712–725. (17) Ryckaert, J.; Ciccotti, G.; Berendsen, J. C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327–341. (18) Nam, K.; Gao, J.; York, D. An efficient linear-scaling Ewald method for long-range electrostatic interactions in combined QM/MM calculations. J. Chem. Theory Comput. 2005, 1, 2–13. (19) Bao, J.; et al. Computational study of bindings of olive leaf extract (OLE) to HIV-1 fusion protein gp41. FEBS Lett. 2007, 581, 2737– 2742. (20) Li, L.; Uversky, V. N.; Dunker, A. K.; Meroueh, S. O. A computational investigation of allostery in the catabolite activator protein. J. Am. Chem. Soc. 2007, 129, 15668–15676. (21) Stoica, I.; Sadiq, S. K.; Coveney, P. V. Rapid and Accurate Prediction of Binding Free Energies for Saquinavir-Bound HIV-1 Proteases. J. Am. Chem. Soc. 2008, 130, 2639–2648. (22) Sanner, M. F.; Olson, A. J.; Spehner, J. C. Reduced surface: an efficient way to compute molecular surfaces. Biopolymers 1996, 38, 305– 320. (23) Sitkoff, D.; Sharp, K. A.; Honig, B. Accurate calculation of hydration free energies using macroscopic solvent models. J. Phys. Chem. 1994, 98, 1978–1988. (24) Wall, L. D.; et al. Binding Constants of Neuraminidase Inhibitors: An Investigation of the Linear Interaction Energy Method. J. Med. Chem. 1999, 42, 5142–5152.

12964

J. Phys. Chem. B, Vol. 114, No. 40, 2010

(25) Cheatham, T. E., III; Srinivasan, J.; Case, D. A.; Kollman, P. A. Molecular dynamics and continuum solvent studies of the stability of polyGpolyC and polyA-polyT DNA duplexes in solution. J. Biomol. Struct. Dyn. 1998, 16, 265–280. (26) Kuhn, B.; Kollman, P. A. Binding of a diverse set of ligands to avidin and streptavidin: an accurate quantitative prediction of their relative affinities by a combination of molecular mechanics and continuum solvent models. J. Med. Chem. 2000, 43, 3786–3791. (27) Cheng, L. S.; et al. Ensemble-Based Virtual Screening Reveals Potential Novel Antiviral Compounds for Avian Influenza Neuraminidase. J. Med. Chem. 2008, 51, 3878–3894.

Wang and Zhang (28) Streicher, H.; Busseb, H. Building a successful structural motif into sialylmimeticsscyclohexenephosphonate monoesters as pseudo-sialosides with promising inhibitory properties. Bioorg. Med. Chem. 2006, 14, 1047– 1057. (29) Wang, J.; Morin, M.; Wany, W.; Kollman, P. A. Use of MMPBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA. J. Am. Chem. Soc. 2001, 123, 5221–5230.

JP1030224