|Home | About | Journals | Submit | Contact Us | Français|
Pandemic (H1N1) influenza poses an imminent threat. Nations have stockpiled inhibitors of the influenza protein neuraminidase in hopes of protecting their citizens, but drug-resistant strains have already emerged, and novel therapeutics are urgently needed. In the current work, the computer program AutoGrow is used to generate novel predicted neuraminidase inhibitors. Given the great flexibility of the neuraminidase active site, protein dynamics are also incorporated into the computer-aided drug-design process. Several potential inhibitors are identified that are predicted to bind neuraminidase better than currently approved drugs.
Influenza is caused by RNA viruses of the family Orthomyxoviridae. While not generally life threatening in healthy adults, the virus occasionally mutates into more deadly forms and has been responsible for several pandemics in the last century. Recently, a new strain of pandemic influenza (H1N1) capable of infecting humans has been identified (Dawood, et al., 2009), with a U.S. hospitalization rate of about 9%. Additionally, a distinct strain of avian influenza (H5N1) arose in 1997 that may cause a similar global pandemic in the future (Abdel-Ghafar, et al., 2008).
In preparation for pandemic influenza, many nations have stockpiled inhibitors of the influenza protein neuaminidase (Oxford, et al., 2004). Following formation, influenza viral particles remain bound to cell membranes via sialic-acid residues. Neuraminidase cleaves these residues, releasing the virus and enabling viral propagation (De Clercq and Neyts, 2007). Neuraminidase is the target of several FDA-approved drugs, including zanamivir and oseltamivir (Oxford, et al., 2004), because it is essential for viral propagation and has a well-conserved active site (Kobasa, et al., 1999). Unfortunately, drug-resistant strains have recently emerged (Kiso, et al., 2004; Beigel, et al., 2005; de Jong, et al., 2005; De Clercq, 2006), and the need for novel inhibitors is great.
Motivated by the urgent need for new influenza therapeutics, we used AutoGrow (Durrant, et al., 2009), a recently developed computer-aided drug-design program, to guide the design of several potential neuraminidase inhibitors predicted to bind better than currently approved drugs.
To account for protein flexibility, we drew upon a molecular dynamics simulation of neuraminidase that has been described previously (Cheng, et al., 2008). Protein conformations extracted from this 40-ns simulation were clustered into 27 groups by root-mean-square-deviation (RMSD) conformational clustering using the gromos clustering algorithm, as implemented in the GROMOS++ analysis software (Daura, et al., 1999; Christen, et al., 2005). In brief, an RMSD distance was calculated for each pair of protein conformations extracted from the MD simulation. Those pairs with associated RMSD distances greater than 1.3 Å were discarded. The single conformation most frequently present in the remaining pairs, together with the other corresponding conformation of each pair, were merged into a list of conformations called the first cluster. The conformations of the first cluster were subsequently removed from the pool of conformations extracted from the MD simulation, and the process was repeated until no conformations remained. The centroid of each cluster was selected, producing an ensemble of 27 unique protein structures representative of the many protein conformations sampled during the simulation.
AutoGrow was run three times, once using a neuraminidase crystal structure (PDB ID: 2HU4, Russell, et al., 2006) as the template protein and twice using each of the top two ensemble structures (Cheng, et al., 2008). In each of these three runs, AutoGrow ran for eight generations, adding fragments to a core scaffold similar to oseltamivir. Each generation initially contained fifty ligands. For each generation after the first, ten primary individuals were taken from the previous generation, based on both the score of the most populated docking cluster and successful active-site binding. An additional twenty “crossovers” and twenty “mutants” were created from these ten primary individuals, subject to the requirement that all compounds contain fewer than one hundred atoms. The first generation initially contained only the scaffold and 49 “mutants,” as no previous generation existed from which “parents” could be drawn for crossover production.
To determine fitness, all AutoGrow-generated ligands were docked into their respective neuraminidase structures (the crystal structure or the two ensemble conformations) using AutoDock 4.0.1 (Morris, et al., 1998), a docking program with a physics-based scoring function that performs well relative to the scoring functions of other similar programs (e.g. DOCK, FleX, and GOLD, Rarey, et al., 1996; Jones, et al., 1997; Ewing, et al., 2001; Bursulaya, et al., 2003). Docking parameters were optimized for the positive-control docking of oseltamivir into the group-1 neuraminidase (N1) crystal structure. The initial AutoDock population size was set to 200 individuals, the maximum number of energy evaluations to 7 × 106, the number of runs to 25, and the RMSD tolerance to 2.0. All other AutoDock parameters were set to the default values. The AutoDock-predicted binding energy was taken to be the energy associated with the most populated AutoDock cluster. AutoDock grids were calculated for regularly spaced points at intervals of 0.375 Å contained within a cube 24.00 Å × 27.00 Å × 24.75 Å, centered on the neuraminidase active site.
To generate novel compounds, AutoGrow drew upon a new fragment library containing 37 637 redundant fragments derived from FDA-approved compounds using an in house script called Fragmentizer. To create this novel fragment library, we first obtained the names of hundreds of FDA-approved compounds by searching Drugs@FDA, supplemented with a list provided by the laboratory of Maurizio Pellecchia. The PDB structures of these compounds were downloaded from drugbank.ca (Wishart, et al., 2006) and filtered to remove those with molecular weight greater than 700 g/mol. After additional processing, 1 174 drugs remained.
For each compound, Fragmentizer first identified all single bonds that could be broken without altering the electronic or geometric configuration of neighboring atoms. The program next generated a second list of all possible bond combinations. Each compound was then decomposed by simultaneously cutting all the bonds of each combination and adding hydrogen atoms to the resulting fragments as needed. Following compound decomposition, all fragments with mass greater than 150 daltons were removed, leaving 37 637 fragments. Redundant fragments were not eliminated. Both Fragmentizer and the novel fragment library derived from FDA-approved compounds can be downloaded from http://www.nbcr.net/software/downloads/virtual_lib/.
As a beta version of AutoGrow was used to generate the ligands, the compounds had to be further processed to correct occasional structural errors. The top ten ligands from each of the three AutoGrow runs were visually inspected. Where the atoms of two distinct fragments were very close, those fragments were bound together to form rings. Where two fragments were mistakenly added via the same scaffold linker hydrogen, extra atoms were removed as needed. Additionally, some sulfur atoms were bound to too many hydrogen atoms. These were eliminated or replaced with oxygen atoms as necessary. Following corrections, each ligand underwent five hundred steps of Cartesian minimization in ICM (Molsoft), a molecular modeling and docking program, prior to being evaluated for drug likeness (Tables 1 and S1).
The relaxed complex scheme (Amaro, et al., 2008a) was used to rescore predicted inhibitors. All compounds were docked into the 27 ensemble configurations using AutoDock 4.0.1 (Morris, et al., 1998). Additionally, six positive controls were included: sialic acid, the natural neuraminidase substrate, in the boat, chair, and twist conformations; zanamivir and oseltamivir, FDA-approved neuraminidase inhibitors; and peramivir, a compound currently in clinical trials. An ensemble-average AutoDock score was calculated for each ligand by averaging the AutoDock scores and weighting according to the cluster population size (Tables 1 and S1):
where Ē is the weighted ensemble-average score, wi is the size of cluster i, and Ei is the AutoDock score of the ligand docked into the centroid of cluster i.
The docking parameters were unchanged, except that the maximum number of energy evaluations was decreased to 5 × 106, and the number of runs was increased to 100.
The top predicted ligand from each of the three AutoGrow runs (compounds 1, 2, and 3; Tables 1 and S1), as judged by the ensemble-average score, was selected for further examination. Each of these three ligands was loaded into ICM (Molsoft), together with the corresponding ensemble member that gave the best score, and redocked using the ICM docking program (Molsoft). This initial docking established a baseline score that was subsequently used to judge whether modifications to the compound improved or reduced binding affinity.
Each of the three ligands then underwent a series of manual modifications, producing fourteen novel drug-like compounds. Modifications were made such that the ICM docking score was maintained, drug-like metrics (molecular weight, number of hydrogen bond acceptors and donors, predicted LogP, etc.) were improved, chirality was reduced, and molecular rigidity was increased. Furthermore, visual inspection led to the addition of some novel groups that AutoGrow had not suggested. These “modified compounds” were rescored using the relaxed complex scheme, with the same AutoDock and AutoGrid parameters used previously.
The influenza virus has caused several pandemics in the last century; recently, a new pandemic strain (H1N1) has been identified (Dawood, et al., 2009). In preparation for pandemic influenza, many nations have stockpiled inhibitors of the influenza protein neuaminidase (Oxford, et al., 2004). Unfortunately, drug-resistant strains have emerged (Kiso, et al., 2004; Beigel, et al., 2005; de Jong, et al., 2005; De Clercq, 2006), and the need for novel neuraminidase inhibitors is great.
Neuraminidase has highly flexible loops adjacent to its sialic-acid binding site. In the first crystal structures of group-1 neuraminidase (N1), the so-called 150-loop adopted an open conformation not seen in previous crystal structures of group-2 proteins (N2) (Russell, et al., 2006). Under different crystallographic conditions, however, the N1 150-loop adopted a closed conformation similar to that of N2 (Varghese, et al., 1983; Baker, et al., 1987). The flexibility of this loop in N1, as well as implications for drug design, have been further characterized via molecular-dynamics simulations (Amaro, et al., 2007; Cheng, et al., 2008).
The great flexibility of the N1 active site defies traditional computer-aided drug-design efforts, which typically focus on static crystal structures and at best account for only limited protein flexibility. To aid future drug-design efforts, Cheng et al. recently performed a 40-ns molecular dynamics simulation of the N1 holo enzyme. 1.6 × 104 protein conformations were extracted at regular intervals and clustered into 27 groups using root-mean-square-deviation (RMSD) conformational clustering. The set of the corresponding 27 distinct centroids, representative of all conformations sampled, is said to constitute an ensemble.
In some fragment-based drug-design strategies, weakly binding molecular fragments, identified experimentally via x-ray crystallography or NMR, or computationally via computer docking, are linked to generate potent composite inhibitors (Rees, et al., 2004). An alternate computational fragment-based growing strategy, exemplified by programs like LUDI (Bohm, 1992; Bohm, 1993, 1994) and AutoGrow (Durrant, et al., 2009), adds small molecular fragments to initial scaffolds known to bind the target protein, with the goal of improving binding affinity. Because fragment-based strategies are combinatorial, a far greater diversity of compounds can be synthesized and tested than would be possible with traditional high-throughput assays.
To generate predicted ligands, we used a beta version of the computer program AutoGrow (Durrant, et al., 2009), an evolutionary algorithm that automates fragment addition to core scaffolds. Each member of a population of AutoGrow-generated compounds was evaluated for binding using AutoDock; the best predicted binders became the founding members of the next generation, wherein fragments were again added/modified. Generation after generation, ligands eventually evolved that were well suited to the template protein structures specified. AutoGrow was initially run against an N1 crystal structure (PDB ID: 2HU4) (Russell, et al., 2006). To account for neuraminidase flexibility, the program was run twice more against the top two representative ensemble conformations.
AutoGrow drew upon a fragment library derived from FDA-approved compounds. The library was generated by decomposing approved drugs into their constituent fragments using an in-house script. As AutoGrow selects fragments from the library at random, redundant fragments were maintained to bias the program towards those fragments most commonly found in FDA-approved drugs. For example, the most common fragment (methane) constituted 11.71% of the library. The ten most abundant fragments of this new library are listed in Table 2. Thirty compounds, ten from each AutoDock run, were considered candidates for further analysis.
To further account for N1 flexibility, the compounds were reranked using the relaxed complex scheme (Amaro, et al., 2008a). In a traditional virtual screen, a docking program is used to predict the ligand energy of binding to a static structure, usually obtained from x-ray crystallography or NMR. The relaxed complex scheme builds upon this traditional methodology by docking candidate ligands into multiple protein conformations extracted from a molecular dynamics simulation. Compounds are then ranked by the ensemble-average predicted binding energy, rather than by the score associated with a single static structure alone. The relaxed complex scheme has been used to identify inhibitors of FKBP (Lin, et al., 2002), HIV integrase (Schames, et al., 2004), and T. brucei RNA editing ligase 1 (Amaro, et al., 2008b).
In addition to rescoring the thirty AutoGrow-generated compounds, we also scored six positive controls: sialic acid, the natural neuraminidase substrate, in the boat, chair, and twist conformations; zanamivir and oseltamivir, FDA-approved neuraminidase inhibitors; and peramivir, a compound currently in clinical trials (Tables 1 and S1). The measured IC50 values of oseltamivir, peramivir, and zanamivir are 0.33 ± 0.27, 0.37 ± 0.26, and 0.57 ± 0.46 nM, respectively (Malaisree, et al., 2008). These IC50 values are roughly three orders of magnitude smaller than the ensemble-average AutoDock-predicted inhibition constants of 2.25, 5.01, and 8.58 μM, respectively. However, there is a clear correlation between prediction and measurement (R2 = 0.92), and, importantly, the relaxed complex scheme ranked the known neuraminidase inhibitors correctly.
The experimentally measured Ki of sialic acid suggests that binding to the natural substrate is much weaker. The measured Ki of 50 μM (Colman, 1994) corresponds well with the ensemble-average AutoDock-predicted inhibition constants of 33.67, 68.02, and 375.48 μM for sialic acid in the boat, twist, and chair conformation, respectively. When the lowest predicted sialic-acid inhibition constant is considered, the correlation between prediction and measurement is further strengthened (R2 = 0.97).
This strong correlation demonstrates that, while the ensemble-average predicted binding energies and inhibition constants are not numerically accurate, the relaxed complex scheme can be used to accurately rank candidate neuraminidase inhibitors. Notably, all thirty AutoGrow-generated compounds ranked better than currently approved neuraminidase inhibitors when the relaxed complex scheme was used; the score of the best novel compound was −12.83 kcal/mol, compared to −8.24 kcal/mol for the best FDA-approved inhibitor, oseltamivir.
Despite the use of fragments derived from FDA-approved drugs, the AutoGrow-generated compounds were not particularly drug like. Consequently, the top ligand from each of the three AutoGrow runs (compounds 1, 2, and 3; Tables 1 and S1), as judged by the ensemble-average score, was selected for further examination. Fourteen drug-like compounds were subsequently generated by making manual modifications to the three lead ligands. Manual modifications were guided by six rules: satisfy Lipinski’s rule of five (Lipinski, et al., 2001), conserve the predicted binding affinity, remove buried unpaired hydrogen bond donors and acceptors, increase molecular rigidity, and reduce chirality.
The AutoDock-generated compounds were generally too large, leading to multiple Lipinski violations (Lipinski, et al., 2001). Lipinski’s rule of five, developed in 2001, states that compounds generally have poor absorption or permeation if they possess two or more of the following properties: more than five hydrogen-bond donors, more than ten hydrogen-bond acceptors, a molecular weight greater than 500 daltons, or a CLogP greater than five. As the AutoGrow-generated compounds were generally large, most contained too many hydrogen bond donors and acceptors and too large a molecular weight. Manual modifications reduced the size of the candidate ligands so they would better conform to Lipinski’s rule of five.
The reduction in size was generally accompanied by a mild drop in the predicted binding affinity. Some interacting groups had to be removed in order to reduce the molecular weight. Additionally, most steps (i.e. chemical modifications) away from the AutoGrow-generated ligand, which had been “optimized” for binding energy, naturally reduced binding affinity. Nevertheless, the compounds were redocked with each manual modification, and modifications that caused precipitous drops in the predicted binding energy were rejected.
Some of the docked AutoGrow-generated compounds had hydrogen-bond donors and acceptors that were buried but unpaired (i.e. “unsatisfied”). The hydrogen-bond donors and acceptors of a fully hydrated ligand are typically satisfied via interactions with the water solvent. A large energy penalty occurs upon binding if unsatisfied hydrogen-bond donors and acceptors are positioned in solvent-inaccessible regions along the protein-ligand interface because a hydrogen bond with the water solvent is lost without the compensatory creation of a protein-ligand hydrogen bond. Some of the manual modifications made served to remove these unpaired hydrogen-bond donors and acceptors.
Additional manual modifications were made to increase compound rigidity. In the unbound state, non-rigid, “floppy” compounds undergo many different conformational transitions as they twist about their rotatable bonds. Upon binding, however, the ligand is stabilized, resulting in a large entropic penalty of binding. For rigid compounds, the difference in entropy between the bound and unbound state is not as great, and so the penalty is less.
These modified compounds were rescored using the relaxed complex scheme (Tables 3 and S2). Ten of the fourteen modified compounds ranked better than FDA-approved neuraminidase inhibitors; the ensemble-average score of the best modified compound was −11.05 kcal/mol, compared to −8.24 kcal/mol for oseltamivir.
The top four modified compounds are all similar, in part because they are built on the same cyclohex-1-enecarboxylate scaffold. While similar to oseltamivir, the novel compounds are decorated with different molecular fragments that, in theory, serve to enhance the binding affinity.
All four compounds have a carboxylate group in the one position, similar to the known neuraminidase inhibitors oseltamivir, paramivir, and zanamivir. An examination of an oseltamivir-neuraminidase crystal structure (Russell, et al., 2006) reveals that this carboxylate group participates in hydrogen bonds with R371 and Y347 (distance cutoff, 3.0; angle cutoff, 30°). In addition, an analysis of the top four modified compounds docked into the 27 ensemble conformations suggests that the carboxylate group also forms hydrogen bonds with R292 and, on rare occasions, with R118 (Figure 1).
All four modified compounds have an (aminomethyl)amino (aminal) group at the five position instead of an amino group, as in oseltamivir. The crystal structure demonstrates that the oseltamivir amine forms hydrogen bonds with D151 and E119. In contrast, an analysis of the top four modified compounds suggests that binding between the aminal fragment and residues D151 and E119 occurs only occasionally. Rather, the distal amine binds instead with E227 and, less frequently, with E277 and the backbone carbonyl of W178. In the case of compounds 5 and 7, the proximal amine is predicted to bind occasionally with E227 as well. Crystal structures indicate that the analogous guanidine groups of peramivir (Zhang, et al., 1992) and zanamivir (Xu, et al., 2008) likewise form hydrogen bonds with the backbone carbonyl of W178 (Figure 1).
Unfortunately, aminal groups are unstable and subject to hydrolysis. Initial efforts to find alternative groups failed; when the aminal was changed to a urea, the ensemble-average predicted binding energy was significantly reduced (data not shown). Another solution may be to replace the aminal with a guanidino group, which is essentially a rigid aminal with an additional amine bound to the bridging carbon atom. We again note that zanamivir and peramivir, neuraminidase inhibitors with IC50’s comparable to that of oseltamivir (Malaisree, et al., 2008), have guanidino groups at the analogous location. Alternatively, the proximal aminal nitrogen atom could be replaced with a carbon atom, though such a change may make synthesis more challenging.
The top four modified compounds contain a 2-methylpropyl group at the four position instead of the acetamido group characteristic of oseltamivir, peramivir, and zanamivir. The crystal structures of oseltamivir, peramivir, and zanamivir bound to neurainidase (Russell, et al., 2006; Xu, et al., 2008) suggest that a hydrogen bond may form between the acetamido carbonyl oxygen atom and R152, though in the case of oseltamivir the hydrogen bond is strained (D-H-A angle of 54°) (Russell, et al., 2006). Additionally, the acetamido methyl group may contribute to the overall binding affinity via hydrophobic interactions with W178. We note, however, that the acetamido amine hydrogen, a hydrogen-bond donor, is buried and unpaired. Like the acetamido group, the 2-methylpropyl group of the top four modified compounds is likewise predicted to form hydrophobic interactions with W178, but without the need for a buried but unpaired hydrogen-bond donor (Figure 1). We note, however, that, in addition to facilitating a hydrogen bond between the carbonyl oxygen atom and the receptor, the amide linker may also simplify chemical synthesis, and so may be necessary on those grounds.
Compounds 4, 5, and 6 have (3-ethyl-2-oxopentan-3-yl)oxy; 2-ethylbutyl; and 2,2-diethyl-3-oxobutyl groups, respectively, at the three position, instead of the pentan-3-yloxy group of oseltamivir. An analysis of the binding poses of the top four modified compounds suggests that the bridging oxygen atom of the oseltamivir pentan-3-yloxy fragment does not contribute significantly to the energy of binding. The crystal structure of oseltamivir bound to neuraminidase reveals no hydrogen bonds with this oxygen atom (Russell, et al., 2006); additionally, an ensemble-wide analysis of compounds 4 and 7, which include the bridging oxygen atom, likewise revealed no hydrogen-bond formation (Figure 1).
The importance of the bridging oxygen atom can be further assessed by comparing the predicted binding energies of compounds that differ only at this location. For example, the ensemble-average predicted binding energies of compounds 5 and 7 (Tables 3), which differ only in the presence or absence of the bridging oxygen atom, are nearly identical. A comparison between compounds 4 and 6 gives a similar result (Tables 3). Though the bridging oxygen is likely energetically unnecessary, we again note that it may greatly facilitate chemical synthesis.
The top four compounds all mimic oseltamivir at the three position in that they contain multiple aliphatic chains. However, compounds 4 and 6 also contain a carbonyl group, a potential hydrogen-bond acceptor. An analysis of these compounds docked into the 27 ensemble conformations revealed that the carbonyl oxygen atom is predicted to form occasional hydrogen bonds with R152. Additionally, the carbonyl oxygen atom of compound 4 forms occasional hydrogen bonds with R292, and the carbonyl oxygen atom of compound 6 forms occasional hydrogen bonds with N294 (Figure 1).
Pandemic influenza (H1N1) poses an imminent threat, and drug-resistant strains have already emerged (Kiso, et al., 2004; Beigel, et al., 2005; de Jong, et al., 2005; De Clercq, 2006). In the current work, the computer program AutoGrow (Durrant, et al., 2009) was used to generate novel predicted neuraminidase inhibitors. Given the great flexibility of the N1 active site, protein dynamics were also incorporated into the computer-aided drug-design process. Several potential inhibitors were identified that are predicted to bind neuraminidase better than current FDA-approved drugs.
The neuraminidase inhibitors oseltamivir and zanamivir, both FDA approved, have similar binding poses (Russell, et al., 2006; Xu, et al., 2008). Despite these similarities, these two compounds have very different resistance profiles. As of 2005, no virus resistant to zanamivir had been isolated from an immunocompetent patient, but resistance to oseltamivir was then emergent (Moscona, 2005). The principal difference between oseltamivir and zanamivir, and the difference likely responsible for their disparate resistance profiles, is the set of molecular fragments used to decorate the central six-member ring. We are therefore hopeful that the novel neuraminidase inhibitors suggested here, with their unique decorating fragments, may likewise serve as scaffolds for future neuraminidase inhibitors against which resistance has not yet developed.
JDD is funded by a Pharmacology Training Grant through the UCSD School of Medicine. This work was carried out with funding from NIH GM31749, NSF MCB-0506593, and MCA93S013 to JAM. Additional support from the Howard Hughes Medical Institute, the National Center for Supercomputing Applications, the San Diego Supercomputer Center, the W.M. Keck Foundation, Accelrys, Inc., the National Biomedical Computational Resource, and the Center for Theoretical Biological Physics is gratefully acknowledged. We would also like to thank Dr. Rommie E. Amaro for helpful discussions and support.
Supporting Information Available