Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Chem Biol Drug Des. Author manuscript; available in PMC 2010 November 16.
Published in final edited form as:
PMCID: PMC2982673

Conformational Dynamics of the Flexible Catalytic Loop in Mycobacterium tuberculosis 1-Deoxy-D-xylulose 5-Phosphate Reductoisomerase


In mycobacteria, the biosynthesis of the precursors to the essential isoprenoids, isopentenyl diphosphate and dimethylallyl pyrophosphate is carried out by the methylerythritol phosphate (MEP) pathway. This route of synthesis is absent in humans, who utilize the alternative mevalonate acid (MVA) route, thus making the enzymes of the MEP pathway of chemotherapeutic interest. One such identified target is the second enzyme of the pathway, 1-Deoxy-D-xylulose 5-phosphate reductoisomerase (DXR). Only limited information is currently available concerning the catalytic mechanism and structural dynamics of this enzyme, and only recently has a crystal structure of Mycobacterium tuberculosis species of this enzyme been resolved including all factors required for binding. Here, the dynamics of the enzyme is studied in complex with NADPH, Mn2+, in the presence and absence of the fosmidomycin inhibitor using conventional molecular dynamics and an enhanced sampling technique, Reversible Digitally Filtered Molecular Dynamics. The simulations reveal significant differences in the conformational dynamics of the vital catalytic loop between the inhibitor-free and inhibitor-bound enzyme complexes and highlight the contributions of conserved residues in this region. The substantial fluctuations observed suggest that DXR may be a promising target for computer-aided drug discovery through the relaxed complex method.


Tuberculosis (TB) is a serious infectious disease caused by the Mycobacterium tuberculosis bacterium. It is a major cause of illness and death with 1.7 million deaths reported globally in 2006(1), and owing to a rise in HIV cases, the neglect of TB control programs and an increase in drug-resistance, the disease has resurged in recent years in well-developed countries and has exacerbated the TB problem in the lesser developed countries(2). Therefore, there is an urgent need for the development of new drugs and suitable therapeutic targets. The enzymes of the methylerythritol phosphate (MEP) pathway (also known as the DOXP or non-mevalonate pathway) are responsible for the biosynthesis of the precursors to isoprenoids (isopentenyl diphosphate (IPP) and its isomer, dimethylallyl (DMAPP)) and have been identified as suitable targets for drug development owing to their vital and diverse functions (e.g. respiration, electron transport, hormone-based signaling and membrane stability(3, 4)) and absence in humans(5, 6).

In most eubacteria, including the causal agents for diseases such as malaria, leprosy and tuberculosis, the MEP pathway is the only route present for the synthesis of isoprenoids (7-10). In archaebateria, fungi and animals the synthesis is exclusively carried out using the mevalonate acid (MVA) pathway (7-9, 11, 12) and in the case of plants, both the MEP and MVA pathways are used, although each is localized in the plastids and cytoplasm respectively (13-17). The existence of these different pathways for the synthesis of IPP and DMAPP between bacteria and humans makes the enzymes of the MEP pathway of particular interest as targets for therapeutic development.

The MEP pathway comprises eight enzymes (18, 19) with 1-deoxy-d-xylulose 5-phosphate reductoisomerase (DXR) being the most studied to date. DXR is involved in the second stage of the pathway, mediating the reversible intramolecular rearrangement and reduction of 1-deoxy-d-xylulose 5-phosphate (DXP) to 2-C-methyl-D-erythritol 4-phosphate (MEP) in the presence of a divalent metal ion (for which Mn2+ has shown to be the most effective (20)) and NADPH.

Drugs such as fosmidomycin, whose structure is similar to the natural substrate, and its analogues have been developed and shown to be efficacious against the E. coli (21) and P. falciparum (22, 23) DXR enzymes. However, these inhibitors are ineffective against the M. tuber strain of the enzyme since mycobacteria are naturally resistant to the majority of antibiotics and chemotherapeutic agents, owing to factors such as the complex nature of the cell wall and the presence of resistance genes (24, 25).

At present, 21 crystal structures of the DXR enzyme are present in the Protein Data Bank (PDB) (26), 11 from E. coli (EcDXR), 8 from M. tuber (MtDXR) and 2 from Z. mobilis (ZmDXR). The 11 E. coli structures have been resolved in the apo-form and in a variety of different complexes involving combinations of the catalytically important divalent metal ion, NADPH cofactor, fosmidomycin inhibitor and bisphosphonate and/or sulfate ions bound in the active site in the place of the natural substrate or inhibitor (27-32). The two ZmDXR enzymes present in the PDB have been resolved in complex with an acetate ion in one structure and with NADPH in the other (33).

Only recently have X-ray crystal structures of MtDXR become available (34, 35) and currently 8 are present in the Protein Data Bank, resolved using a new cloning, expression and purification technique, making the M. tuber enzyme more receptive to structural analysis (34). These eight crystal structures are of the wild-type and double mutant (D151N and E222Q) forms of the enzyme in the apo-state and in complex with combinations of Mn2+, NADPH, fosmidomycin and SO42- bound.

In each of the crystal structures from M. tuber, the enzyme is observed to form an asymmetric homodimer (figure 1). Each 413-residue monomer comprises three domains, an N-terminal domain (binds NADPH), a catalytic domain (harbors active site) and C-terminal domain (figure 1). The N- and C-terminal domains are positioned in a V-shape with the catalytic domain located at the cleft and are shown to share a 40 and 35 % sequence identity with the DXR enzymes of E. coli and Z. mobilis. The catalytic domain shows higher conservation of its residues with a 50 % similarity of sequence identity compared with the E. coli and Z. mobilis (34).

Figure 1
Crystal structure of asymmetric homodimer of MtDXR (PDB ID: 2JCZ). Monomer A has been coloured according to the following structural features: blue: N-terminal, red: C-terminal, orange: catalytic domain, yellow: extended loop of catalytic domain, purple: ...

The active site comprises three regions into which the substrate binds. The phosphate moiety of the substrate binds into a positively charged site, a hydrophobic pocket interacts with the substrate backbone, and an amphipathic region binds the hydroxamate portion (30). A cluster of conserved acidic residues, Asp151, Glu153 and Glu222 gives rise to the binding site of the catalytically vital divalent metal ion. The octahedral binding geometry of this metal ion in MtDXR is revealed in a crystal structure of the enzyme in complex with NADPH, Mn2+ and fosmidomycin, resolved by Henriksson et al. (34)

Near the active site lies a flexible loop structure (residues 198-209 in MtDXR), which, in conjunction with the C-terminal and catalytic domain are indicated by crystal structure conformations to undergo conformational change on binding of a substrate or inhibitor (27). This loop has been resolved in two major conformations, open and closed, and variations in between. As with the EcDXR structures (30, 32), where the M. tuber enzyme is bound to a substrate (or fosmidomycin) and NADPH cofactor, the loop is observed in a closed “lid” conformation, reaching over the active site with loop residues interacting with the substrate (or fosmidomycin) (34). The NADPH cofactor has been found to be essential in the tight binding of the substrate or inhibitor owing to the contribution of the nicotinamide ring to the formation of the hydrophobic binding pocket (30). In DXR crystal structures where a substrate or fosmidomycin is bound but NADPH is missing, the loop is shown to be open (29) or in one crystal structure containing a SO42- ion in the place of fosmidomycin and no NADPH, an intermediate conformation is observed (35). In a study by Henriksson et al, they note that the loop of the EcDXR structures appears to be able to attain a more closed conformation compared with that of MtDXR and associate this with either the species variation or an artifact of the different crystallization conditions used (34).

The purpose of the loop closure is to provide a solvent-shielded cavity for substrate processing which additionally provides suitable dielectric properties for the binding and processing of a small and highly charged substrate (30). The rearrangement of the loop to this closed conformation provides an intramolecular hydrophobic surface of approximately 80 % (30).

Several highly conserved residues in the loop and active site are proposed to play an important role in binding, such as loop residues, His200 (His209 in EcDXR), Trp203 (Trp212), Met205 (Met214) and Met267 (Met276) (30). The histidine is found to form a hydrogen bond with the phosphate group of the substrate or fosmidomycin inhibitor, and is thought to be responsible for fixing the loop in the closed conformation over the active site (28). Mutagenesis studies of H209Q have confirmed the importance of this residue, finding a significant decrease in catalytic activity of the mutant EcDXR enzyme (36). Trp203, Met205 and Met267 are further conserved residues which assume orientations which shield the active site from the surrounding bulk solvent on substrate and inhibitor binding (37). Further to this, the β-indole portion of Trp203 lies parallel and provides a key interaction with the backbone of the substrate. A mutation study by Fernandes et al. also suggests the residue to play a role in the discrimination of substrates after it was found that the tryptophan residue caused a significant decrease in catalytic activity of the enzyme with substrate analogues (38).

The most complete picture of the interactions in the active site of MtDXR is provided by the recent Henriksson et al. structure (34) (PDB ID: 2JCZ) possessing the NADPH cofactor, Mn2+ ion and fosmidomycin. The structure is of the homodimer, however only monomer A includes all the factors required for binding and shows the catalytic loop in a closed conformation, whereas monomer B is fixed in an open conformation due to crystal contacts.

As mentioned previously, experimental work suggests an induced fit mechanism of substrate/fosmidomycin binding involving significant conformational change of a loop located near the active site although the details of the mechanism are not clear. In this study, the enhanced sampling technique, reversible digitally filtered molecular dynamics (RDFMD) simulations, using the recently resolved quaternary MtDXR crystal structure (PDB ID: 2JCZ) are used to describe the dynamics of the enzyme, focusing on the crucial catalytic loop region. In addition, simulations using conventional MD have been used to supplement and compare the sampling with RDFMD simulations. The results describe significant differences between the fosmidomycin bound and unbound states, providing an insight to events occurring on binding.


The starting structure (PDB ID: 2JCZ) was taken from the Protein Data Bank. This crystal structure is of the asymmetric homodimer, resolved with monomer A possessing an Mn2+ ion, NADPH and fosmidomycin bound in the active site. Monomer B is not so well resolved and the density of a large portion of NADPH is poor and only a sulfate is found bound in the binding site in the place of fosmidomycin. As explained in the introduction section, monomer B is thought to be significantly affected by crystal contacts and the active site is held in the open conformation. In this study, simulations were carried out on the asymmetric homodimer but owing to the reasons just described, results are only reported on monomer A.

This structure has been resolved to 2.05 Å as a truncated form of the enzyme with the C-terminus missing 24 residues. However, this truncation does not affect the catalytic activity and the enzyme is found to be fully functional in this state (34). Additionally, the first 10 residues are absent from the structure owing to poor electron density.

Simulations have been carried out with the MtDXR crystal structure in the presence (referred to here as DXR-FMN, the same naming given in the study detailing the crystallization of this structure (34)) and absence of the fosmidomycin (DXR-MN) inhibitor in monomer A. In both cases, NADPH and Mn2+ are present. The DXR-MN structure was generated by deleting the fosmidomycin inhibitor from the crystal structure.

Simulations were performed within the active pH range of the enzyme, reported to be between pH 7.4 and 7.9 (39). Standard protonation states were set for all ionizable residues with the exception of Glu and Asp residues lying in close vicinity to the Mn2+ ion. These residues were changed to their unprotonated form. The WHATIF program (40) was used to add polar hydrogens and to check the DXR-FMN and DXR-MN structures. The AMBER utility XLEAP (41) was used to add other hydrogen atoms and to solvate the systems with a minimum distance of 12 Å from the protein, in a box of 26047 TIP3P water molecules. Thirty-one and twenty-eight Na+ counterions were added to neutralize the overall charge of the DXR-FMN and DXR-MN systems respectively.

All simulations, unless otherwise stated have been carried out using the NAMD molecular dynamics package (42) and the AMBER99 forcefield (43). The parameters for NADPH (44) and Mn2+ (45) were obtained from the AMBER parameter database (

Minimization was carried out in stages, starting with the protein only (5,000 steps), followed by solvent (30,000 steps), ions (1,000 steps), solvent and ions (20,000 steps), NADPH (2500 steps), protein, ions and NADPH (3000) and finally the entire system (90,000 steps), giving a total of 101,500 steps. The two default minimization algorithms present in the NAMD package were used, the line search and the conjugate gradient algorithms.

The minimized system was heated to 300 K in the NVT ensemble. The procedure employed a Langevin thermostat with a 10 ps-1 damping parameter and a 2 fs timestep. The heating was carried out gradually in stages at 50 K intervals, each interval being 20,000 steps long. Following this, a further 50,000 steps in the NVT ensemble was carried out using a 5 ps-1 thermostat damping parameter to control the temperature of the system.

Equilibration simulations, using a Nose-Hoover barostat in the NPT ensemble were then carried out for 50,000 steps, with a target pressure of 1 atm. A decay parameter of 100 fs and a piston period of 200 fs were used. A further 50,000 steps were run, with a decay parameter of 300 fs and a piston period of 500 fs. For each system, two equilibration runs were carried out to produce two starting structures with different velocities. The final equilibrated systems had box dimensions of approximately x: 84 y: 86 z: 123 Å.

The solvation of the binding site cavity was checked and found to be adequate using the minimization and equilibration protocol described here. Proper solvation of cavities is critical, as vacuums can form and artificial artifacts in the dynamics of the system can arise.

For each of the two systems, two MD production simulations of 20 ns in length have been carried out. All production MD simulations were run in the NVT ensemble using a 2 fs timestep, a Langevin thermostat with a 5 ps-1 damping parameter at a temperature of 300 K. Periodic boundary conditions were used along with a Particle Mesh Ewald treatment of electrostatic interactions, using an interpolation order of 6 and switching function applied to the Lennard-Jones interactions between 9 Å and the 10.5 Å cutoff. PME gridsizes of 90 × 90 × 135 Å were used, similar values to those of the boxsizes. SHAKE (46) was applied to all bonds containing hydrogen, using a tolerance of 10-8 Å.

RDFMD simulation details

Reversible Digitally Filtered Molecular Dynamics (RDFMD) enhances conformational change through amplification of the low frequency motions of specific structural regions of a protein. Prior to this study, the method has been successfully applied to a number of protein systems, including E. coli dihydrofolate reductase (47), HIV-1 integrase (48) and HIV-1 protease (49). The details of the method have been described previously in the literature and a protocol of parameters for use on regions of interest in proteins has also been developed previously (47, 50).

The parameters are heavily interrelated and some are system dependent and therefore, a suitable set of parameters has been optimized through trial and error for the study of the MtDXR systems simulated here. A range of protocols are used so to investigate the reproducibility of results when using this non-equilibrium method. These include the use of a digital filter designed to amplify frequencies between 0-100 cm-1 using 201 coefficients and amplification factors of 1.6, 1.8, 2.0, 2.2 and 2.4. A temperature cap of 800 K and a delay between filter applications of either 50 or 100 steps were also employed. The filter sequences were separated by 4 ps of molecular dynamics simulation in the NVT ensemble. This is sufficient time for the system temperature to return to 300 K and it is during this period of time that conformations for analysis are generated. The final results are taken from piecing together 100 individual 4 ps MD runs, totaling 400 ps for each RDFMD simulation.

It is important to note that RDFMD is a non-equilibrium method. For the purposes of this study, it is not designed to generate a thermodynamically representative ensemble; instead, it is a method for generating perturbed structures which would suggest possible binding geometries that would be thermodynamically evaluated in a subsequent step.

Since the dynamics of the catalytic loop is thought to play a fundamental role in the activity of the enzyme, the residues of this region (residues 198-209) have been selected as the target region of the filter in the RDFMD simulations.

Simulations were run in the NVT ensemble using the Langevin thermostat with a 5 ps-1 damping parameter.

A total of 12 RDFMD simulations have been carried out using the same two starting structures for each system, as described for MD, and using a range of different filter delay and amplification factors listed here.

Cluster analysis

Clustering was used to study the different loop conformations sampled. The RMSD-based GROMOS algorithm (51) as implemented in g_cluster was used and analysis was carried out on two trajectories, one formed from the concatenation of MD (structures extracted at 10 ps intervals) and RDFMD trajectories of DXR-FMN. The same was carried out for the simulations of the DXR-MN to generate the second trajectory. For each pair of structures, a least-square translational and rotational fit was performed using all the Cα atoms of monomer A, and the RMSD for the Cα atoms of the loop residues were calculated. A cutoff of 1.2 Å was chosen after evaluation of the cluster populations at a range of cutoff values. All structures within this cutoff were assigned as neighbors with the structure possessing the largest number of neighbors being assigned as the centre of the cluster. This structure was removed from the pool along with its neighbors and the process repeated until all the structures are assigned to clusters.

Principal component analysis

Principal component analysis (PCA) is a method that simplifies a dataset by reducing its dimensions to the representative principal components hopefully without losing any of the important characteristics. It captures the variance of the dataset in terms of principal components, with the first principal component (eigenvector with the highest associated eigenvalue) contributing the most to the variance in the dataset. The calculation and analysis of the principal components were carried out by the g_covar and g_anaeig modules available in the GROMACS 3.3 (52).

The analysis was performed by initially removing rotational and translational motion from the conformations and building a covariance matrix of the positional fluctuations obtained from the MD and RDFMD trajectories for Cα atoms for the residues of monomer A. Diagonalization of the matrix generates a set of eigenvectors and associated eigenvalues, which define a new set of generalized coordinates, and the eigenvalues are sorted in order of decreasing value. It is often found that the first few eigenvectors describe the significant motions of the protein with further eigenvectors often corresponding to smaller random fluctuations. In this work, the inner-products of eigenvector sets are calculated to assess the similarity between the eigenvectors.


Two 20 ns MD simulations (referred to as MD1 and MD2) and twelve 0.4 ns RDFMD simulations (referred to as R1-12) were carried out for both the DXR-FMN and DXR-MN systems. As mentioned in the methods section, simulations were carried out on the enzyme as a dimer but results reported here focus on monomer A. This monomer is significantly less affected by crystal contacts and has improved resolution especially in the loop region and for the NADPH cofactor where poor electron density is observed for monomer B.

As a check of the stability, the RMSD of the Cα atoms of the protein from its equilibrated structure was determined for each of the simulations (figure 2), as well as monitoring of the secondary structure and potential energy (data not shown), all of which confirmed protein stability throughout the simulations.

Figure 2
RMSD of Cα atoms from the equilibrated starting structure over the length of MD1 simulations of DXR–MN (red line) and for DXR–FMN (black line). Results of the other simulations were observed to be similar and are not shown here ...

For the RDFMD simulations, the peptide bonds of the catalytic loop residues, the target residues of the digital filter, were additionally monitored. Owing to an improbable trans to cis isomerization of one peptide bond of the catalytic loop residues, one of the RDFMD trajectories of the DXR-FMN binary complex was discarded, resulting in a total of eleven RDFMD simulations for this system.

Global Motion

Analysis of the Root Mean Squared Fluctuations (RMSF) of the Cα atoms of each of the residues of monomer A reveals the catalytic domain, particularly in the region of the catalytic loop, to be one of the most flexible regions of the enzyme in simulations of both the DXR-FMN and DXR-MN complexes (figure 3). This result is in agreement with experimental studies, which found an important feature of this enzyme to be its intrinsic flexibility (27, 28). The enzyme is proposed to undergo a large conformational change on binding to a substrate or inhibitor, which includes an increase in order of the loop from its more disordered state in the absence of the substrate or inhibitor.

Figure 3
Average RMSF of Cα atoms of monomer A calculated from (A) 12 RDFMD simulations of DXR–MN, (B) 2 MD simulations of DXR–MN (C) 11 RDFMD simulations of DXR–FMN and (D) 2 MD simulations of DXR–FMN.

The variability in RMSF between the simulations for each system, especially between the two MD simulations of the DXR-FMN complex (figure 3 (d)) makes it difficult to conclude that these are differences in flexibility between the inhibitor-bound and inhibitor-free complexes. However, comparison of the 12 RDFMD simulations for the two complexes does indicate reduced flexibility of the catalytic loop where fosmidomycin is bound.

Note that no comparisons between the RMSF of MD and RDFMD simulations have been made. The RDFMD simulations clearly indicate an increase in flexibility, especially in the region of the catalytic loop, but this region is the target of the digital filter, so an increase in flexibility is expected.

Conformations of the catalytic loop sampled by simulations of the DXR-FMN and DXR-MN complexes

To summarize the range of catalytic loop conformations sampled by MD and RDFMD simulations of the DXR-FMN and DXR-MN complexes, cluster analysis was employed, using snapshots generated by the simulations (as described in the methods section).

Cluster analysis of DXR-MN structures

Clustering of the concatenated MD and RDFMD trajectories using an RMSD-based cutoff of 1.2 Å resulted in a total of 45 clusters from RDFMD simulations and 5 clusters from MD simulations. Figure 4 displays the sampling and the transitions between the clusters as a function of time for the simulations, and figure 5 shows representative loop conformations taken from each of the more significant clusters, which together encompass ~90 % of the total conformations sampled by the MD and RDFMD simulations.

Figure 4
Sampling of clusters over 40 ns MD (MD1 and MD2) and 4.8 ns RDFMD (R1–12) simulations.
Figure 5
Representative conformations of the catalytic loop from the clustering of snapshots generated from 4.8 ns of RDFMD and 40 ns of MD simulation. Numbers in brackets depict the % of the total number of snapshots in each cluster. Mn2+ is shown in green van ...

The greater number of clusters identified as a result of RDFMD simulations indicates that increased conformational sampling was achieved using this methodology compared with conventional MD. Approximately 98 % of the snapshots from 40 ns of MD simulation were represented in the first cluster (figure 5 C1), where the loop is covering the active site in a closed conformation, close to that of the starting conformation of the loop. The use of the enhanced sampling technique, RDFMD resulted in the sampling of a variety of conformations with eight clusters incorporating >80 % of the total number of snapshots generated by the simulations. The clusters show the loop in several different conformations, described here as closed, intermediate and open. Closed structures are defined here as when the loop is highly ordered and covering the binding site. Open structures reveal access to the active site, and intermediate structures show a small lifting motion of the loop but the active site remains covered. The most closed-type structures are observed in clusters 1 and 6 (figure 5: C1 and C6) with structures of cluster 1, where ~50 % of the snapshots are clustered being the most similar to that of the crystal structure starting conformation. The most open structures are found in clusters 5 and 8, where the loop is shown to lift in an upwards direction on one side of the loop and to twist, folding away from the active site, both motions resulting in exposure of the binding cavity. Other conformations observed are intermediate between closed and open (clusters 2, 3, 4 and 7), where the loop is shown to have slightly lifted on one side in the region of the His200 residue and is not as ordered as the compact closed loop, but still remains covering the active site. One such conformation of MtDXR has been resolved using X-ray crystallography (PDB ID:2C82). Henriksson et al. resolved the enzyme with a sulfate ion in the place of the substrate in the binding site with no NADPH cofactor or divalent metal ion bound (35). This loop conformation shows similarity with the intermediate conformations sampled by RDFMD simulations, especially with the representative structure of cluster 2. However, the relevance of this crystal structure has been questioned as the structure was crystallized at a pH where the enzyme is not active and therefore, the protonation states and interactions of some residues will differ at the active pH (30).

Cluster analysis of DXR-FMN structures

The same clustering analysis carried out using snapshots generated from MD and RDFMD simulations of the DXR-FMN quaternary complex demonstrated significantly less mobility of the catalytic loop. As seen for the clustering results of DXR-MN, RDFMD simulations of DXR-FMN sampled an increased number of loop conformations (23 clusters) compared with MD simulations (5 clusters) (figure 6), although the total number of clusters (23 clusters) found is reduced compared with simulations where the inhibitor was absent (45 clusters) and ~90 % of the total number of generated snapshots are represented by the first four clusters (compared with 8 when fosmidomycin was absent) (figure 7). Figure 6 shows the conformations sampled by MD1 and MD2 to lie in different clusters (cluster 1 and cluster 2 respectively), both of which show the loop to cover the active site with the loop of cluster 2 being slightly lifted in the vicinity of His200, whereas structures of cluster 1 are more similar to the starting closed crystal structure. This loop conformation sampled predominately throughout MD2 is only briefly sampled in RDFMD simulations, in simulation R5. RDFMD simulations of DXR-FMN are shown to have mainly sampled closed conformations and are members of cluster 1. The significantly smaller sized third cluster, sampled only by RDFMD simulations shows an intermediate conformation with loop covering the active site with the tip of the loop moved upwards.

Figure 6
Sampling of clusters over 40 ns MD (MD1 and MD2) and 4.4 ns RDFMD (R1–11) simulations of DXR–FMN.
Figure 7
Representative conformations of the catalytic loop from the clustering of snapshots generated from 4.4 ns of RDFMD and 40 ns of MD simulation. Numbers in brackets depict the % of the total number of snapshots in each cluster. Mn2+ is shown in green van ...

In summary, clustering analysis shows that neither the MD nor RDFMD simulations see the opening of the loop as observed from the analysis of simulations of the DXR-MN ternary complex.

Importance of conserved loop residues His200 and Trp203

As discussed in the introduction section, His200 and Trp203 are two conserved loop residues, which have been highlighted in experimental studies to play important roles in the catalytic function of the DXR enzyme and their mutation to other residues has been found to significantly decrease the catalytic activity of the enzyme.

Overall, evaluation of the structures from clustering of the snapshots generated by simulations of DXR-MN and DXR-FMN agrees with experimental results. In the absence of fosmidomycin, the Trp203 and His200 residues show an increase in flexibility and can be found far from the active site owing to the absence of interaction with the inhibitor, whereas, in the most closed structures, both residues are positioned close to the active site in orientations agreeing with their proposed roles. His200 forms a hydrogen bond with the phosphonate group of fosmidomycin and the β-indole ring of Trp203 stacks onto the backbone of the inhibitor.

Examination of the occurrence of the intermediate conformations in the simulations indicates this intermediate conformation to occur through a lifting motion of the loop in the vicinity of the His200 residue with the loop remaining over the active site. At these times, the His200 is frequently found away from the active site, even where fosmidomycin is bound, breaking the interaction with the phosphonate group seen to be present in the more closed structures, whereas the Trp203 residue remains close to the backbone of the inhibitor. The crystal structure of MtDXR (PDB ID: 2C82) resolved in an intermediate conformation when bound to a sulfate ion also sees the His200 residue positioned away from the active site whilst the Trp203 residue is positioned towards the active site.

Overall, the results appear to suggest the hydrogen bond between the phosphonate group of the inhibitor and His200 is required for a fully closed structure but the loop can still remain covering the active site in an intermediate conformation even in the absence of this bond, although it is not known whether this is an active conformation. At these times, as seen in the simulations of DXR-MN, the catalytic loop is prevented from fully opening and folding back away from the active site by other loop residues, including Trp203 which maintains its interaction in the presence of the inhibitor in the active site.

The observations of the interactions present between these two residues and the fosmidomycin inhibitor agrees with conclusions reached by experimental mutagenesis studies. Mutation of Trp203 in EcDXR, to residues possessing aliphatic side-chains resulted in a significant decrease in the specificity and catalytic rate of the enzyme, confirming the importance of the interaction between the aromatic side-chain of Trp and the backbone of the substrate (38). The same interaction between Trp203 and the fosmidomycin, a structural analogue of the natural substrate, was observed in the simulations presented in this study.

Similarly, a mutagenesis study of His209Gln in EcDXR also observes a decrease in catalytic rate and is thought to, in part, to play a role in binding the substrate (36).

In the absence of these interactions with the substrate, it may be that the loop is unable to close adequately to form the required environment for substrate processing. However, with the unavailability of crystal structures containing these mutations, it cannot be confirmed.

Opening and closing of the catalytic loop

To measure the extent of opening occurring in simulations, the RMSD of the loop Cα atoms of the MD and RDFMD trajectories from the crystal structure conformation of the loop in the open (monomer B in PDB ID: 2JCZ) and closed (monomer A PDB ID: 2JCZ) states were calculated.

Figure 8 shows the RMSD as a function of time over the 40 ns MD simulation of DXR-FMN and DXR-MN, and for both systems there is no significant deviation in RMSD and no notable conformational change occurs, as also shown by the cluster analysis. The application of digital filters to the loop residues in the RDFMD simulations of the DXR-MN complex indicates several significant opening events (figure 9) and the conformations of the loop where the RMSD fluctuations are greatest are overlaid in figure 10. Although similarity is observed between the open loop conformations sampled by the RDFMD simulations and that of the open crystal structure (PDB ID: 2JCZ), where the loop is folded back, revealing the active site to the bulk solvent, none of the RDFMD simulations sample the open loop to such an extent as seen in the crystal structure. This maybe due to the effect of crystal contacts stabilizing the open conformation of monomer B in the crystal structure, as mentioned previously.

Figure 8
RMSD of Cα atoms of the catalytic loop against the crystal structure (PDB ID: 2JCZ). Top two lines in each plot show the RMSD against the open loop (monomer B of 2JCZ) and bottom two lines show the RMSD against the closed loop (monomer A of 2JCZ) ...
Figure 9
RMSD of catalytic loop residues measured against crystal structure of open loop (bold) and closed loop (dotted) over the 12 and 11 RDMFD simulations of (A) DXR–MN (B) DXR–FMN respectively.
Figure 10
Loop conformations of three opening events observed in RDFMD simulations of DXR–MN (see numbered labels in Figure 9A). Coloring of loop structures (corresponding label in Figure 9A shown in brackets): red (1), orange (2), yellow (3). Open and ...

Smaller fluctuations in the RMSD are seen for RDFMD simulations of DXR-FMN (figure 9(b)), none of which lead to significant opening and the loop remains covering the active site in all simulations in closed or intermediate conformations.

Principal component analysis

To identify correlated motions that are significant to loop opening, PCA was used (as described in the methods section) to study the RDFMD simulations where the catalytic loop is shown to open the greatest. For each of these simulations, PCA was carried out using the Cα atom trajectories of all residues of monomer A and figure 11 displays a representative matrix generated from the analysis of the inner-products of the simulations. A value of unity represents identical eigenvectors whereas a value of zero is given for orthogonal eigenvectors.

Figure 11
Representative inner-product matrix of the first 10 eigenvectors for two RDFMD simulations of DXR–MN that observe opening of the catalytic loop. Similarity between eigenvectors demonstrated by shading. High to low similarity represented by dark ...

The first eigenvector contributes ~50 % of the total motion in each of these simulations and inner-product analysis shows high similarity between the first eigenvectors of the simulations, with an average inner-product of 0.468 (SD 0.117), thus showing the first eigenvector to dominate the total motion of this enzyme.

Visualization of the first eigenvector (figure 12) is shown using a porcupine plot where each cone is placed at the location of a Cα atom and represents the direction (length of cone) and relative magnitude of motion, scaled to the number of atoms in the analysis (53). This eigenvector captures the lifting of the loop in the region of the His200 residue, which results in the folding back of the loop. In conjunction with this loop opening, the N- and C-terminal domains are observed to move in opposite directions resulting in the enlargement of the active site area and causing NADPH to move away from the binding site in the more open structure. The increased distance between the N- and C-terminals has also been noted in experimental studies which have used the distance as a rough guide to the extent of opening seen (34). Movement of the extended loop and rest of the catalytic domain are also observed. Further eigenvectors contribute significantly less to the total motion of the subunit and display considerably reduced similarity, as shown in Figure 11, owing mostly to their representation of smaller localized motions of the loop and the rest of the enzyme.

Figure 12
View from the top of monomer A of MtDXR. The direction and length of the cone pointing from each of the Cα of the average DXR structure show the direction and magnitude of motion captured by the first eigenvector. N-terminal highlighted in red, ...


MD and RDFMD simulations have been used to investigate the conformational dynamics of the MtDXR enzyme in the presence and absence of the fosmidomycin inhibitor, both complexes possessing the essential Mn2+ ion and NADPH cofactor.

MD simulations sampled mainly the closed conformation for both complexes and no significant conformational change was observed.

The use of the enhanced sampling technique, RDFMD greatly enhanced the conformational sampling of the enzyme and the simulations observed a significantly larger range of conformations compared with simulations using conventional MD, especially in the case of the ternary complex, where the inhibitor is absent.

RDFMD simulations of the fosmidomycin-bound complex show the catalytic loop of the enzyme to predominantly sample closed-type conformations where the loop covers the active site. Two highly conserved residues, Trp203 and His200 are found close to the active site in orientations in agreement with experimental studies, with the β-indole ring of Trp203 interacting with the backbone of fosmidomycin and the existence of a hydrogen bond between His200 and the phosphonate group of the inhibitor.

In the absence of fosmidomycin, RDFMD simulations show the enzyme to display increased flexibility, sampling a range of conformations of the closed, intermediate and open enzyme. Where opening occurs, PCA and clustering of the simulations reveal the loop to open through a lifting and twisting motion, causing the loop to fold back away from the active site, exposing the area to the bulk solvent. At these times, the loop residues, Trp203 and His200 are frequently found far from the active site.

Clustering shows a significant proportion of conformations sampled by the inhibitor-free simulations, and a small proportion of inhibitor-bound simulations, to be intermediate between open and closed, where the loop is slightly lifted in the vicinity of the His200 residue but still remains covering the active site. Examination of these conformations finds the His200 residue is able to move away from the active site, even where fosmidomycin is bound, breaking the interaction with the phosphonate group of the inhibitor. In simulations of the inhibitor bound MtDXR simulations, the Trp203 residue maintains its interaction with the backbone of fosmidomycin and together with other residue interactions, keeps the loop over the active site. A similar loop conformation is seen in the crystal structure of MtDXR bound with a sulfate ion, however, the relevance of this structure has been questioned and it is not known whether this conformation represents an active form of the enzyme and has not been seen when bound to the natural substrate. The functioning enzyme may require the full loop closure, therefore confirming the necessity of the hydrogen bond between His200 and the substrate. The intermediate structure may represent the first stages of the opening of the loop.

PCA has revealed that, on loop opening, the N- and C-terminals move in opposing directions, opening up the active site area with simultaneous motion of the extended loop of the catalytic domain. The results are in agreement with experimental evidence, which reports the purpose of the domain motion and loop closure is to form a tight and specific binding site. The rearrangement of the loop to the closed conformation provides a solvent-shielded cavity, increasing the hydrophobic nature of the active site. This may alter the relationship between charged residues within the binding site, providing a microscopic environment with suitable dielectric conditions for the binding and processing of the small and highly-charged substrate

An improved understanding of the dynamics of the catalytic loop and the interactions that occur could prove useful in the discovery of new inhibitors. The loop is postulated to be crucial in the successful binding of inhibitors and substrates. The simulations presented here observe the closed loop to provide a small and hydrophobic binding site, which would enable the formation of specific interactions required for substrate processing. A more open loop is coupled with a more solvent-exposed binding site. The RDFMD simulations here have sampled a range of conformations, some of which have not been reported previously, and have provided an insight into the role of the highly conserved His200 and Trp203 loop residues and the interconversion from closed to more open conformations via an intermediate conformation. The variety of conformations sampled suggests the MtDXR enzyme maybe a suitable target for computer-aided drug discovery techniques such as the relaxed complex method (54-56), which has previously been used to account for receptor flexibility in systems such as HIV-1 integrase (55) and avian influenza neuraminidase (57). Such techniques examine numerous configurations of the target, enabling the incorporation of the intrinsic flexibility of the binding site and increasing the likelihood of the discovery of novel and chemically diverse inhibitors.


This work was supported in part by the NSF, NIH, HHMI, CTBP and NBCR.


1. WHO. Global tuberculosis control: surveillance, planning, financing: WHO report 2008. World Health Organization; 2008. Global tuberculosis control: surveillance, planning, financing: WHO report 2008.
2. Ebrahim GJ. Drug resistance in tuberculosis. J Trop Pediatr. 2007;53:147–9. [PubMed]
3. Sacchettini JC, Poulter CD. Biochemistry - Creating isoprenoid diversity. Science. 1997;277:1788–9. [PubMed]
4. Lange BM, Wildung MR, McCaskill D, Croteau R. A family of transketolases that directs isoprenoid biosynthesis via a mevalonate-independent pathway. Proc Natl Acad Sci U S A. 1998;95:2100–4. [PubMed]
5. Swanson KM, Hohl RJ. Anti-cancer therapy: Targeting the mevalonate pathway. Curr Cancer Drug Targets. 2006;6:15–37. [PubMed]
6. Jomaa H, Wiesner J, Sanderbrand S, Altincicek B, Weidemeyer C, Hintz M, Turbachova I, Eberl M, Zeidler J, Lichtenthaler HK, Soldati D, Beck E. Inhibitors of the nonmevalonate pathway of isoprenoid biosynthesis as antimalarial drugs. Science. 1999;285:1573–6. [PubMed]
7. Boucher Y, Doolittle WF. The role of lateral gene transfer in the evolution of isoprenoid biosynthesis pathways. Mol Microbiol. 2000;37:703–16. [PubMed]
8. Kuzuyama T. Mevalonate and nonmevalonate pathways for the biosynthesis of isoprene units. Biosci Biotechnol Biochem. 2002;66:1619–27. [PubMed]
9. Lange BM, Rujan T, Martin W, Croteau R. Isoprenoid biosynthesis: The evolution of two ancient and distinct pathways across genomes. Proc Natl Acad Sci U S A. 2000;97:13172–7. [PubMed]
10. Rohdich F, Kis K, Bacher A, Eisenreich W. The non-mevalonate pathway of isoprenoids: genes, enzymes and intermediates. Curr Opin Chem Biol. 2001;5:535–40. [PubMed]
11. Beytia ED, Porter JW. Biochemistry of polyisoprenoid biosynthesis. Annu Rev Biochem. 1976;45:113–42. [PubMed]
12. Edwards PA, Ericsson J. Sterols and isoprenoids: Signaling molecules derived from the cholesterol biosynthetic pathway. Annu Rev Biochem. 1999;68:157–85. [PubMed]
13. Lichtenthaler HK. The 1-deoxy-D-xylulose-5-phosphate pathway of isoprenoid biosynthesis in plants. Annual Rev Plant Physiol Plant Mol Biol. 1999;50:47–65. [PubMed]
14. Lichtenthaler HK, Rohmer M, Schwender J. Two independent biochemical pathways for isopentenyl diphosphate and isoprenoid biosynthesis in higher plants. Physiol Plant. 1997;101:643–52.
15. Eisenreich W, Schwarz M, Cartayrade A, Arigoni D, Zenk MH, Bacher A. The deoxyxylulose phosphate pathway of terpenoid biosynthesis in plants and microorganisms. Chem Biol. 1998;5:R221–R33. [PubMed]
16. Eisenreich W, Rohdich F, Bacher A. Deoxyxylulose phosphate pathway to terpenoids. Trends Plant Sci. 2001;6:78–84. [PubMed]
17. Rohmer M. The discovery of a mevalonate-independent pathway for isoprenoid biosynthesis in bacteria, algae and higher plants. Nat Prod Rep. 1999;16:565–74. [PubMed]
18. Rohdich F, Wungsintaweekul J, Fellermeier M, Sagner S, Herz S, Kis K, Eisenreich W, Bacher A, Zenk MH. Cytidine 5′-triphosphate-dependent biosynthesis of isoprenoids: YgbP protein of Escherichia coli catalyzes the formation of 4-diphosphocytidyl-2-C-methylerythritol. Proc Natl Acad Sci U S A. 1999;96:11758–63. [PubMed]
19. Hunter WN. The non-mevalonate pathway of isoprenoid precursor biosynthesis. J Biol Chem. 2007;282:21573–7. [PubMed]
20. Takahashi S, Kuzuyama T, Watanabe H, Seto H. A 1-deoxy-D-xylulose 5-phosphate reductoisomerase catalyzing the formation of 2-C-methyl-D-elythritol 4-phosphate in an alternative nonmevalonate pathway for terpenoid biosynthesis. Proc Natl Acad Sci U S A. 1998;95:9879–84. [PubMed]
21. Kuzuyama T, Shimizu T, Takahashi S, Seto H. Fosmidomycin, a specific inhibitor of 1-deoxy-D-xylulose 5-phosphate reductoisomerase in the nonmevalonate pathway for terpenoid biosynthesis. Tetrahedron Lett. 1998;39:7913–6.
22. Missinou MA, Borrmann S, Schindler A, Issifou S, Adegnika AA, Matsiegui PB, Binder R, Lell B, Wiesner J, Baranek T, Jomaa H, Kremsner PG. Fosmidomycin for malaria. Lancet. 2002;360:1941–2. [PubMed]
23. Lell B, Ruangweerayut R, Wiesner J, Missinou MA, Schindler A, Baranek T, Hintz M, Hutchinson D, Jomaa H, Kremsner PG. Fosmidomycin, a novel chemotherapeutic agent for malaria. Antimicrob Agents Chemother. 2003;47:735–8. [PMC free article] [PubMed]
24. Jarlier V, Nikaido H. Mycobacterial cell wall - structure and role in natural-resistance to antibiotics. FEMS Microbiol Lett. 1994;123:11–8. [PubMed]
25. Nguyen L, Thompson CJ. Foundations of antibiotic resistance in bacterial physiology: the mycobacterial paradigm. Trends Microbiol. 2006;14:304–12. [PubMed]
26. Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, Shindyalov IN, Bourne PE. The Protein Data Bank. Nucleic Acids Res. 2000;28:235–42. [PMC free article] [PubMed]
27. Reuter K, Sanderbrand S, Jomaa H, Wiesner J, Steinbrecher I, Beck E, Hintz M, Klebe G, Stubbs MT. Crystal structure of 1-deoxy-D-xylulose-5-phosphate reductoisomerase, a crucial enzyme in the non-mevalonate pathway of isoprenoid biosynthesis. J Biol Chem. 2002;277:5378–84. [PubMed]
28. Yajima S, Nonaka T, Kuzuyama T, Seto H, Ohsawa K. Crystal structure of 1-deoxy-D-xylulose 5-phosphate reductoisomerase complexed with cofactors: Implications of a flexible loop movement upon substrate binding. J Biochem. 2002;131:313–7. [PubMed]
29. Steinbacher S, Kaiser J, Eisenreich W, Huber R, Bacher A, Rohdich F. Structural basis of fosmidomycin action revealed by the complex with 2-C-methyl-D-erythritol 4-phosphate synthase (IspC) - Implications for the catalytic mechanism and anti-malaria drug development. J Biol Chem. 2003;278:18401–7. [PubMed]
30. Mac Sweeney A, Lange R, Fernandes RPM, Schulz H, Dale GE, Douangamath A, Proteau PJ, Oefner C. The crystal structure of E. coli 1-deoxy-D-xylulose-5-phosphate reductoisomerase in a ternary complex with the antimalarial compound fosmidomycin and NADPH reveals a tight-binding closed enzyme conformation. J Mol Biol. 2005;345:115–27. [PubMed]
31. Yajima S, Hara K, Sanders JM, Yin FL, Ohsawa K, Wiesner J, Jomaa H, Oldfield E. Crystallographic structures of two bisphosphonate : 1-deoxyxylulose-5-phosphate reductoisomerase complexes. J Am Chem Soc. 2004;126:10824–5. [PubMed]
32. Yajima S, Hara K, Iino D, Sasaki Y, Kuzuyama T, Ohsawa K, Seto H. Structure of 1-deoxy-D-xylulose 5-phosphate reductoisomerase in a quaternary complex with a magnesium ion, NADPH and the antimalarial drug fosmidomycin. Acta Crystallogr Sect F Struct Biol Cryst Commun. 2007;63:466–70. [PMC free article] [PubMed]
33. Ricagno S, Grolle S, Bringer-Meyer S, Sahm H, Lindqvist Y, Schneider G. Crystal structure of 1-deoxy-D-xylulose-5-phosphate reductoisomerase from Zymomonas mobilis at 1.9-angstrom resolution. Biochim Biophys Acta. 2004;1698:37–44. [PubMed]
34. Henriksson LM, Unge T, Carlsson J, Aqvist J, Mowbray SL, Jones TA. Structures of Mycobacterium tuberculosis 1-deoxy-D-xylulose-5-phosphate reductoisomerase provide new insights into catalysis. J Biol Chem. 2007;282:19905–16. [PubMed]
35. Henriksson LM, Bjorkelid C, Mowbray SL, Unge T. The 1.9 angstrom resolution structure of Mycobacterium tuberculosis 1-deoxy-D-xylulose 5-phosphate reductoisomerase, a potential drug target. Acta Crystallographr D Biol Crystallogr. 2006;62:807–13. [PubMed]
36. Kuzuyama T, Takahashi S, Takagi M, Seto H. Characterization of 1-Deoxy-D-xylulose 5-phosphate reductoisomerase, an enzyme involved in isopentenyl diphosphate biosynthesis, and identification of its catalytic amino acid residues. J Biol Chem. 2000;275:19928–32. [PubMed]
37. Singh N, Avery MA, McCurdy CR. Toward Mycobacterium tuberculosis DXR inhibitor design: homology modeling and molecular dynamics simulations. J Comput Aided Mol Des. 2007;21:511–22. [PubMed]
38. Fernandes RPM, Phaosiri C, Proteau PJ. Mutation in the flexible loop of 1-deoxy-D-xylulose 5-phosphate reductoisomerase broadens substrate utilization. Arch Biochem Biophys. 2005;444:159–64. [PubMed]
39. Dhiman RK, Schaeffer ML, Bailey AM, Testa CA, Scherman H, Crick DC. 1-Deoxy-D-xylulose 5-phosphate reductoisomerase (IspC) from Mycobacterium tuberculosis: towards understanding mycobacterial resistance to fosmidomycin. J Bacteriol. 2005;187:8395–402. [PMC free article] [PubMed]
40. Vriend G. WHAT IF - A molecular modeling and drug design program. J Mol Graph. 1990;8:52–6. [PubMed]
41. Case DA, Darden TA, Cheatham TE, III, Simmerling CL, Wang J, Duke RE, Luo R, Merz KM, Wang B, Pearlman DA, Crowley M, Brozell S, Tsui V, Gohike H, Mongan J, et al. Amber 8. University of California; San Francisco: 2004.
42. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel RD, Kale L, Schulten K. Scalable molecular dynamics with NAMD. J Comp Chem. 2005;26:1781–802. [PMC free article] [PubMed]
43. Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA. A 2nd generation force-field for the simulation of proteins, nucleic acids, and organic-molecules. J Am Chem Soc. 1995;117:5179–97.
44. Holmberg N, Ryde U, Bulow L. Redesign of the coenzyme specificity in L-lactate dehydrogenase from Bacillus stearothermophilus using site-directed mutagenesis and media engineering. Protein Eng. 1999;12:851–6. [PubMed]
45. Bradbrook GM, Gleichmann T, Harrop SJ, Habash J, Raftery J, Kalb J, Yariv J, Hillier IH, Helliwell JR. X-ray and molecular dynamics studies of concanavalin-A glucoside and mannoside complexes - Relating structure to thermodynamics of binding. J Chem Soc Faraday Trans. 1998;94:1603–11.
46. Ryckaert JP, Ciccotti G, Berendsen HJC. J Comput Phys. 1977;23:327–41.
47. Wiley AP, Swain MT, Phillips SC, Essex JW, Edge CM. Parameterization of reversible digitally filtered molecular dynamics simulations. J Chem Theory Comput. 2005;1:24–35. [PubMed]
48. Williams SL, Essex JW. The study of the conformational dynamics of the catalytic loop of WT and G140A/G149A HIV-1 Integrase core domain using Reversible Digitally Filtered Molecular Dynamics. 2008 unpublished work. [PubMed]
49. Wiley AP, Williams SL, Essex JW. Conformational motions of HIV-1 Protease identified using Reversible Digitally Filtered Molecular Dynamics. 2008 unpublished work. [PubMed]
50. Phillips SC, Swain MT, Wiley AP, Essex JW, Edge CM. Reversible digitally filtered molecular dynamics. J Phys Chem B. 2003;107:2098–110.
51. Daura X, Gademann K, Jaun B, Seebach D, van Gunsteren WF, Mark AE. Peptide folding: When simulation meets experiment. Angew Chem Int Ed Engl. 1999;38:236–40.
52. Van der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJC. GROMACS: Fast, flexible, and free. J Comput Chem. 2005;26:1701–18. [PubMed]
53. Barrett CP, Hall BA, Noble MEM. Dynamite: a simple way to gain insight into protein motions. Acta Crystallogr D Biol Crystallogr. 2004;60:2280–7. [PubMed]
54. Lin JH, Perryman AL, Schames JR, McCammon JA. The relaxed complex method: Accommodating receptor flexibility for drug design with an improved scoring scheme. Biopolymers. 2003;68:47–62. [PubMed]
55. Schames JR, Henchman RH, Siegel JS, Sotriffer CA, Ni HH, McCammon JA. Discovery of a novel binding trench in HIV integrase. J Med Chem. 2004;47:1879–81. [PubMed]
56. Lin JH, Perryman AL, Schames JR, McCammon JA. Computational drug design accommodating receptor flexibility: The relaxed complex scheme. J Am Chem Soc. 2002;124:5632–3. [PubMed]
57. Cheng LS, Amaro RE, Xu D, Li WW, Arzberger PW, McCammon JA. Ensemble-based virtual screening reveals potential novel antiviral compounds for avian influenza neuraminidase. J Med Chem. 2008;51:3878–94. [PMC free article] [PubMed]