Utilizing a Dynamical Description of IspH to Aid in the Development of Novel Antimicrobial Drugs

The nonmevalonate pathway is responsible for isoprenoid production in microbes, including H. pylori, M. tuberculosis and P. falciparum, but is nonexistent in humans, thus providing a desirable route for antibacterial and antimalarial drug discovery. We coordinate a structural study of IspH, a [4Fe-4S] protein responsible for converting HMBPP to IPP and DMAPP in the ultimate step in the nonmevalonate pathway. By performing accelerated molecular dynamics simulations on both substrate-free and HMBPP-bound [Fe4S4]2+ IspH, we elucidate how substrate binding alters the dynamics of the protein. Using principal component analysis, we note that while substrate-free IspH samples various open and closed conformations, the closed conformation observed experimentally for HMBPP-bound IspH is inaccessible in the absence of HMBPP. In contrast, simulations with HMBPP bound are restricted from accessing the open states sampled by the substrate-free simulations. Further investigation of the substrate-free simulations reveals large fluctuations in the HMBPP binding pocket, as well as allosteric pocket openings – both of which are achieved through the hinge motions of the individual domains in IspH. Coupling these findings with solvent mapping and various structural analyses reveals alternative druggable sites that may be exploited in future drug design efforts.


Introduction
In the past couple decades, antimicrobial drug resistance has risen dramatically and greatly hampered the efficacy of currently available therapies for bacterial and malarial infections [1][2][3][4][5][6][7][8][9]. Whereas (multiple-)drug-resistant bacterial infections are a ubiquitous problem, affecting both the Western world and developing nations, the burdens of malaria fall disproportionately on the poorest regions of the world, with over 219 millions cases and 666,000 deaths reported in 2010 [3]. Beyond the common problems associated with decreased lifetimes for drug efficacy due to rapid development of resistance [1,2,5,6,9], advances in the fight against bacterial and malarial infections have also been plagued by diminished attention from major pharmaceutical companies toward the development of new therapies and drugs [4,5,10]. Consequently, there is urgent need for the development of new drugs with novel modes of action, for administration either independently or in combination with established regimen, both to combat bacterial and malarial infections, as well as to address the propensity of each for rapidly developing drug resistance [1,4,6,7,9].
The nonmevalonate pathway for isoprenoid biosynthesis has recently been revealed as a novel target for both antibacterial and antimalarial drugs. Isoprenoids comprise essential metabolites derived from the 5-carbon biomolecules, isopentenyl diphosphate (IPP) and dimethylallyl diphosphate (DMAPP, Figure 1), examples of which include sterols that provide structural support to membranes, chlorophylls used in photosynthesis, and quinones that participate in electron transport chains [11][12][13][14]. In contrast, animals acquire IPP and DMAPP in a distinctive manner via a mevalonate-dependent pathway. Given this metabolic difference, the proteins involved in the nonmevalonate pathway provide novel targets for the development of antibacterial and antimalarial drugs that are both broadly specific to pathogenic species such as H. pylori, M. tuberculosis and P. falciparum and without known human analogs [15][16][17][18].
The ultimate step in the nonmevalonate pathway is the generation of IPP and DMAPP through a 2-electron reductive dehydroxylation of (E)-1-hydroxy-2-methyl-but-2-enyl pyrophosphate (HMBPP) by IspH, a [4Fe-4S] protein ( Figure 1) [11,[19][20][21]. The catalytic mechanism of IspH has been a topic of great debate, largely due to uncertainties introduced by the iron-sulfur cluster [18,22]. Initial structures of IspH from Aquifex aeolicus [23] and Escherichia coli [24] solved by X-ray crystallography resemble cloverleaves and comprise three sequentially different domains with pseudo-C3 symmetry, each tethered to a [Fe 3 IspH) assumes an open conformation, with a 10620 Å cavity where the HMBPP molecule is expected to bind at the cluster [23]. In contrast to the A. aeolicus crystal structure, the [Fe 3 S 4 ] + E. coli counterpart is closed around an inorganic diphosphate molecule (PP i ) that sits in the vicinity of the centrally located [Fe 3 S 4 ] + cluster. Various conserved polar and charged residues, including Glu-126, Thr-167, Asn-227, His-41, His-74, His-124, Ser-225, Ser-226 and Ser-269 (E. coli numbering scheme), coordinate the PP i molecule, likely via hydrogen bonding or salt bridge interactions [24]. The orientations of these conserved residues in the E. coli structure are distinct from their A. aeolicus counterparts due to a tilt of a single domain that enables co-localization of charged and polar residues around the PP i in the case of the former.
While results from electron paramagnetic resonance (EPR) spectroscopy have shown [Fe 3 S 4 ] + IspH to be catalytically active [25], reconstituted IspH displays EPR and Mossbauer signatures of a [Fe 4 S 4 ] 2+ cluster [26,27]. Groll and co-workers provide further support for the catalytically relevant form of IspH containing a [Fe 4 S 4 ] 2+ cluster with their work in crystallizing IspH in the presence of its substrate, HMBPP. This HMBPPbound crystal structure (PDB ID: 3KE8, henceforth referred to as [Fe 4

S 4 ] 2+
(closed, HMBPP-bound) IspH) assumes a closed conformation having a domain tilt similar to that of the [Fe 3 S 4 ] + E. coli structure, with HMBPP bound via its terminal hydroxyl moiety to an unliganded iron of a [Fe 4 S 4 ] 2+ cluster ( Figure 2) [28]. The coordination sphere of the HMBPP ligand is virtually identical to the inorganic diphosphate molecule, while its terminal hydroxyl moiety interacts with Glu-126, Thr-167 (E. coli numbering) and an ordered water molecule to make a hydrogen bond network that is proposed to facilitate proton transfer during catalysis [28]. While these structural data provide a good picture of the [Fe 4 S 4 ] 2+ IspH structure with HMBPP bound, the structure of the 4Fe-form in the absence of substrate, as well as a detailed understanding of how IspH changes conformation upon ligand binding, are not fully understood.
Drawing from insight gained from the aforementioned structural work, as well as various spectroscopic and mutational studies, multiple groups have contributed to drug discovery efforts on the IspH target [29][30][31][32][33][34]. To the best of our knowledge, IspH inhibitor development has fallen under two classes: (1) HMBPP analogues [29][30][31] and (2) pyridine or alkenyl/alkynyl diphosphates and bisphosphonates [32][33][34]. In the case of HMBPP analogues, inhibitor binding emulates the natural substrate, while leveraging improved interactions with the Fe-site (e.g. binding of a thiol instead of an alcohol) [30,31]. Alternatively, Oldfield and coworkers have created novel inhibitors of IspH by utilizing olefinic and pyridine groups to form p/s ''metallacycle'' complexes and g 1 -complexes, respectively, coupling these metal binding groups to phosphate skeletons that preserve the hydrogen bond and salt bridge interactions present in IspH-HMBPP complexation [32][33][34]. These initial drug discovery efforts may be enhanced, both in terms of finding new lead compounds and developing already discovered leads, by obtaining a better description of the IspH binding pocket and possible allosteric sites that may be targeted.
Given that there exists no high-resolution structural data for substrate-free, [Fe 4 S 4 ] 2+ IspH, this work employs accelerated molecular dynamics (aMD) simulations to describe the dominant conformations available to IspH having a fourth iron atom in the absence of HMBPP. Characterization of these dominant conformations reveals an expanded binding pocket and allosteric sites that may be targeted with future rational drug design efforts. Additional attention is directed toward understanding how IspH

Author Summary
Drug resistance has recently entered into media conversations through the lens of MRSA (methicillin-resistant Staphylococcus aureus) infections, but conventional therapies are also failing to address resistance in cases of malaria and other bacterial infections, such as tuberculosis.
To address these problems, we must develop new antibacterial and antimalarial medications. Our research focuses on understanding the structure and dynamics of IspH, an enzyme whose function is necessary for the survival of most bacteria and malaria-causing protozoans. Using computer simulations, we track how the structure of IspH changes in the presence and absence of its natural substrate. By inspecting the pockets that form in the normal motions of IspH, we propose a couple new routes by which new molecules may be developed to disrupt the function of IspH. It is our hope that these structural studies may be precursors to the development of novel therapies that may add to our current arsenal against bacterial and malarial infections.
dynamics are altered upon ligand binding, allowing us to propose a mechanism for how IspH-HMBPP complexation is achieved.

RMSD and visual analyses of aMD simulations of open, substrate-free IspH
Consistent with the nomenclature used by Gräwert, et al. [28], descriptions of IspH from this point forward will use the nomenclature D1, D2 and D3 to describe the domains containing residues 14-96; 97-193; and 194-281, 1-13, respectively (A. aeolicus numbering, Figure 2). We perform 36100 ns aMD simulations of [Fe 4

S 4 ] 2+
(open, substrate-free) IspH, starting from the A. aeolicus crystal structure with a fourth iron modeled into the cluster, as described in the Methods. All trajectories are aligned to the [Fe 3 S 4 ] + (open, substrate-free) IspH crystal structure by the backbone atoms of all D1 residues, since these residues display significantly lower fluctuation throughout the simulation than those in D2 and D3 [23]. The root-mean-square deviation (RMSD) for the backbone atoms of all residues after alignment is given in Figure 3a. From this RMSD analysis, it is apparent that each independent trajectory samples conformational space differently. The large changes in RMSD correspond to opening and closing motions of the D2 and D3 domains, providing a more dynamic description of the [Fe 4

Docking of HMBPP to open IspH
Using Schrodinger's Glide program [35][36][37], we dock HMBPP to the unique iron site in IspH. Docked poses are filtered applying knowledge from experiment that the terminal alkoxide/alcohol group of HMBPP directly chelates the apical Fe site [26,28,38,39]. The docked pose used in our MD studies is found by constraining the position of the terminal alkoxide moiety to within a 2.5 Å radius of the apical iron. While the orientation of the PP i moiety in our docked pose differs from the [Fe 4 S 4 ] 2+ (closed, HMBPP-bound) IspH crystal structure (3KE8) [28], it is worth mentioning that the cyclic structure of HMBPP observed in the crystal structure likely results from ''induced fit'' effects, with polar and charged groups closing around the PP i moiety. Given these effects are absent from our docking procedure, we use the Glide geometry as a starting point for elucidating how open, substrate-free IspH responds to the formation of an encounter complex with HMBPP bound to its unliganded Fe. (open, substrate-free) IspH crystal structure, with the RMSD of all backbone atoms to the crystal structure given in Figure 3b. Both seeds one and three (Figure 3b, black and blue, respectively) approach an RMSD of ,8-10 Å , with respect to the crystal structure. This jump occurs rapidly for seed three (in the first 20 ns of simulation), while seed one only appears to approach this level in the last 10  (closed, HMBPP-bound) IspH crystal structure ( Figure S1).   (closed, HMBPP-bound) IspH crystal structure ( Figure 4). When considering the structures in this closed cluster, it is notable that the ligand does not form a ring structure consistent with its pose in the crystal structure. Nevertheless, the closing of the D2 and D3 domains around the substrate is consistent with the [Fe 4 S 4 ] 2+ (closed, HMBPP-bound) experimental reference [28]. A more detailed inspection of the HMBPP environment in a representative structure from this closed cluster reveals the three key active site histidines, as well as the conserved Thr-165, Thr-166, Glu-126, Ser-221, Asn-223 and Ser-265 forming contacts with HMBPP that appear identical to those seen in the [ crystal structure ( Figure 4b). These findings demonstrate that aMD simulations have effectively captured the closing of loops from D2 and D3 around HMBPP-confirming earlier hypotheses for how conformational change occurs upon substrate binding [28].

RMSD and visual analyses of aMD simulations of apo-IspH with HMBPP docked into the active site
Inconsistent with the RMSD results for seed three, both seeds one and two (black and red, Figure 3b Figure S1). Similar to what is seen in the [Fe 4 S 4 ] 2+ /HMBPP (open, docked) aMD simulations, we note that HMBPP never fully reaches its ring conformation seen crystallographically [28]. From these simulations, it is evident that substrate-bound IspH, once folded around HMBPP, has less conformational space accessible to it and does not access open states. Assessing sampling using principal component analysis

In constructing principal component (PC) space from all [Fe 4 S 4 ] 2+
(open, substrate-free) and [Fe 4 S 4 ] 2+ /HMBPP (open, docked) simulations, as described in the Methods, we observe that the first two principal components account for 83% of the variance. Using Bio3D [40], the motions that correspond to movement along PC1 and PC2 are visualized ( Figure S2) and are shown to correspond to opening and closing motions achieved through the hinge-like properties of the loops that connect D3 to D1 and D2 and D2 to D1 and D3, as suggested by Groll and co-workers [28].
All  (Figure 5b). Three minima are apparent in the projections, one centered near the substrate-free crystal structure and two near the HMBPP-bound crystal structures that differ slightly in the specific contacts made between the protein and ligand. The extent to which bound-HMBPP restricts IspH dynamics is highlighted from the projection of the [Fe 4 S 4 ] 2+ / HMBPP (closed) simulations onto PC space. When simulated from a closed conformation, it is clear that bound-HMBPP effectively locks the protein in a closed conformation, unable to access open states-evident by a well present only around the closed, HMBPPbound IspH crystal structure (Figure 5c).

RMSF analysis shows a decrease in fluctuation accompanies ligand binding
Combining the three trajectories for each individual system simulated, we performed a root-mean-square fluctuation (RMSF) analysis to quantify the extent to which each residue fluctuates in the IspH (white diamond, PDB ID: 3KE8) are also projected onto PC space [23,28]. Numbers in (A) correspond to POVME volumes (Å 3 ), as described in the text. Obtaining an understanding of local phenomena driving conformational change upon HMBPP binding Changes in various peptide dihedral angles (phi, psi and chi) typically accompany global changes in protein conformation [42,43]. In other words, certain dihedral angles may select for specific conformations in proteins [42,43]. Recently, McClendon et al. contributed a method that quantifies differences in probability distributions of protein dihedral angles between a reference and altered state of a protein by using an expansion of the Kullback-Leibler (KL) Divergence [42]. This application assigns a value for the ''mutual divergence'' of each residue-a measure of the extent to which the distributions of dihedral angles differ between the two states. Using the [Fe 4

S 4 ] 2+
(open,substrate-free) simulation as a reference, we compute ''mutual divergence'' values upon substrate binding to IspH using the MutInf suite of programs [42,44] in an attempt to isolate local changes in protein structure that give rise to globally different conformational ensembles between the [Fe 4
A visual representation of ''mutual divergence'' values is provided in Figure 7a, with the highest scoring residues shown in Table 1. We present these data together with measures of sequence conservation, computed as Shannon entropy [40,45,46], as both these metrics are suggested to highlight residues of functional importance [42,45]. The link between sequence conservation and functionality is obvious-residues that are highly conserved are usually conserved for some purpose, e.g. to bestow certain structural features to a protein or to participate in catalysis. Similarly, residues whose conformations change dramatically upon some natural perturbation to the system, ligand binding in our case, are likely responsible for the functionality of that protein.
Consequently, we propose that residues that both are highly conserved and display high ''mutual divergence'' upon ligand binding are critical to the structure and function of IspH. Most other residues with high ''mutual divergence'' can be characterized by one of two distinct environments in the protein: either (a) coordinating HMBPP when it is bound (e.g. His-42, His-124, Asn-223 and Ser-265); or (b) structurally flanking the thiolates that anchor the [Fe 4 S 4 ] 2+ cluster to the protein (as is the case for Phe-12, which is adjacent to Cys-13).
High mutual divergence is seen for residues that occupy the first coordination shell of HMBPP when it is bound. These residues assume different conformations based upon whether they are  (open,substrate-free) state. These differences derive from the reorientation of these residues about the PP i of HMBPP in order to participate in hydrogen bonds or salt bridges.
The other class of residues with high ''mutual divergence'' reside adjacent to the thiolates tethered to the [Fe 4 S 4 ] 2+ cluster. Phe-12 exemplifies this finding, in that it maintains altered w/y angle distributions, contingent on whether HMBPP is bound ( Figure S3 Applying knowledge of the expanded binding pocket from apo simulations to motivate drug discovery efforts

Clustering of the [Fe 4 S 4 ] 2+
(open,substrate-free) IspH aMD simulations reveals dominant structures with substrate pockets of differing volumes and chemical environments. Using representative structures from each of the clusters, we investigate the druggability of these different pockets by performing solvent mapping with the FTMAP program [48].
Taking the fragment positions as they are docked by FTMAP into each representative structure from the [Fe 4

S 4 ] 2+
(open,substrate-free) simulations, we synthesize information regarding where the docked fragments congregate by generating a probe occupancy map for IspH. Probe occupancy is highest at the pocket corresponding to the substrate-binding site (Figure 8a, Figure S4). In the more voluminous clusters, as well as the most dominant cluster, probes expand beyond the HMBPP-binding site at the iron, into a crevice between D1 and D3 (Figure 8a, Figure S4). This finding suggests that inhibitors capable of occupying this expanded pocket while locking the protein in a state that is more open with respect to the [Fe 3

S 4 ] +
(open, substrate-free) IspH crystal structure may provide a feasible route toward novel inhibitor design.
An unanticipated finding from solvent mapping concerns the side of IspH opposite the substrate-binding pocket. When the protein opens fully, as seen in the [Fe 4

S 4 ] 2+
(open,substrate-free) simulations, the hinge-like quality of the interface between D3 and D1/D2 hyperextends, creating a druggable pocket found opposite the side of the HMBPP-binding site (Figure 8b, black rectangle; Figure S4). When the hinge is opened, this pocket occupies a POVME-measured volume of 330-500 Å 3 and accommodates a variety of polar and nonpolar probes. This result, stemming from the opening motions intrinsic to substratefree, [Fe 4 S 4 ] 2+ IspH, may provide an allosteric target for inhibitor design.

Application of the aMD method to sample conformational space in both [Fe 4 S 4 ] 2+
(open,substrate-free) and [Fe 4 S 4 ] 2+ / HMBPP (closed) states of IspH increases our understanding of how HMBPP binding affects IspH structure and dynamics, as well as highlights alternative routes for the design of novel IspH inhibitors.
In regard to IspH dynamics, our [Fe 4 S 4 ] 2+ /HMBPP (open,docked) aMD simulations are able to capture the closing event that accompanies ligand binding in two out of three simulations. In these simulations, residues in D1 that are needed to coordinate HMBPP (His-42 and His-74) are already properly positioned to interact with the pyrophosphate tail of HMBPP, whereas residues in D2 and D3 that coordinate HMBPP require domain motions to bring them in proximity of the substrate. Once D2 and D3 close around HMBPP, it is apparent from our PCA that IspH is unable to reopen, with the fluctuations of residues from D2 and D3 largely suppressed as these domains engage in multiple electrostatic and hydrogen bond interactions with HMBPP (e.g. His-124, and Ser-226). These observations underscore the suggestions by others that both electron addition to the substrate and changes in IspH crystal structure, indicating that substrate binding allows IspH to sample a closed state that is inaccessible in the absence of HMBPP.
In our simulations of the substrate-free state, IspH accesses both open and closed conformations. Since closed states preexist in the substrate-free ensemble, it is tempting to suggest that conformational selection (CS) [49] is responsible for ligand recognition in IspH. Following the logic of Sullivan and Holyoak [50], however, induced fit (IF) likely better describes the conformational changes occurring upon ligand binding since HMBPP cannot actually bind to the closed state that preexists in substrate-free IspH (due to occlusion of the active site by D2 and D3) [51]. We propose that ligand binding may still be described as a combination of CS and IF, where the ligand initially selects open conformations for formation of an encounter complex. Once initially bound, HMBPP induces closure of D2 and D3 via motions that are also intrinsic to IspH in the absence of ligand. This ligand recognition mechanism, drawing from both CS and IF, is not unique to IspH, but rather gives further support to the suggestions of others that ligand binding can contain elements of both CS and IF [52][53][54].
Returning state with the results from our KL divergence analysis and FTMAP solvent mapping of IspH, we can build on the work of others [29][30][31][32][33][34] in suggesting a novel framework for future IspH inhibitor design. HMBPP binding to IspH can be regarded the first step in the catalytic process vital to most microbes for production of IPP and DMAPP. Preventing this binding event is thus the goal of competitive inhibitor development. From our KL divergence analysis, we find that in addition to conserved residues that coordinate HMBPP upon its binding, residues that are adjacent to thiolate residues achieve high ''mutual divergence'' scores due to their distinct dihedral distributions when IspH is open and closed. Given its position adjacent to the fully conserved Cys-193, Asn-194 likely coordinates the hinge motions of D3 that give way to the necessary closing events that accompany HMBPP binding. Preventing the closing of the D3 hinge and, consequently, locking the Asn-194 backbone dihedrals in their disordered, open conformations may provide a novel mode of inhibiting IspH.
From our aMD simulations, two differing mechanisms for disrupting the hinge motions of D3 are apparent. The first targets the outward motion of D3 from the HMBPP binding site that creates an enlarged cavity that extends from the active site to the interface between D3 and D1 (Figure 8a). Either design of larger competitive inhibitors that interact with the apical iron and the D3/D1 interface or design of ligands that interact allosterically with the D3/D1 interface could successfully exploit the enlarged pocket on the active site side of IspH. Alternatively, the presence of an allosteric pocket opposite the side of the HMBPP binding site may be targeted for inhibitor design (Figure 8b, Figure S4). Both these proposed sites for inhibitor design are ''hot spots'' found by solvent probes with FTMAP. Noting that probe occupancy correlates well with sequence conservation as measured by Shannon entropy (r = 0.49) provides further support for these suggested modes of inhibition.
Given the documented difficulties of rational drug design for metalloproteins, notably from a computational perspective [56], allosteric sites that do not require a detailed description of metal binding (e.g. orbital interactions, polarization and charge transfer) are highly desirable if existent. Furthermore, it has been shown that perturbations to allosteric networks in redox-active metalloproteins may affect the redox potential of these proteins and, consequently, alter their activities [57]. These factors motivate us to include the different pockets revealed by aMD simulations, particularly those that may provide routes to noncompetitive inhibition, in future computer aided drug design workflows.

Conclusion
Using aMD simulations, we are able to capture the closing event that accompanies the binding of HMBPP to IspH when starting from the substrate-free crystal structure. Drawing from PCA and visual analyses of the different trajectories considered, we propose that ligand binding occurs via a combination of induced fit and conformational selection. We note that a single dihedral angle, the y angle in Asn-194, selects for either open or closed conformations of IspH, the latter being achieved via a crank motion that draws D3 inward to corral the active site. Furthermore, our aMD simulations reveal both an expanded active site pocket encompassing a crevice between D1 and D3, as well as an allosteric pocket between D1 and D3 on the side opposite the substrate binding pocket that may be utilized for the design of novel IspH inhibitors. Using the Amsterdam Density Functional program [58], a model [Fe 4 S 4 ] 2+ cluster is geometry optimized using broken symmetry density functional theory (BS-DFT) [59,60] at the OLYP/TZP level of theory [61,62]. With the Gaussian 09 suite of programs [63], we optimize the geometry of HMBPP and compute the electrostatic potentials of both geometry optimized HMBPP and the model [Fe 4 S 4 ] 2+ cluster using MK radii [64] at the HF/6-31G(d) level of theory. The antechamber program [65] in the AmberTools 13 suite of programs [66] is then used to apply the restrained electrostatic potential (RESP) procedure to derive point charges for use in MD simulations. In the case of the [Fe 4 S 4 ] 2+ cluster, parameters for nonbonded terms are taken from the AMBER GAFF force field [67], and bonds and angles between atoms are implicitly accounted for by harmonic restraints applied to these terms, using parameters from the [Fe 4 S 4 ] 2+ (closed, HMBPP-bound) IspH crystal structure [28]. For HMBPP, all force field parameters are taken from the AMBER GAFF force field [67]. All charge and nonbonded parameters, as well as, a more detailed discussion of the ligand parameterization process, are included in the Supporting Information.

System preparation for production molecular dynamics simulations
Hydrogens are added using PDB2PQR [68,69], with protonation states assigned using the PROpKa program [70]. In our setup, His-42 and His-124 are set to their imidazolium states, and Glu-126 is protonated. Following hydrogen addition, the protein systems are minimized for 2000 steps in the gas phase using the sander module in AMBER12 [66], to remove problematic steric clashes. The systems are solvated in a box of TIP3P waters [71] that extends 12 Å beyond the closest solute atom, with counterions added to enforce electroneutrality. Non-water bonds to hydrogen atoms are constrained using the SHAKE algorithm [72], while the O-H bonds in water are constrained using the SETTLE algorithm [73]. All protein force field parameters are taken from the AMBER ff99SB force field [74], while the ligand parameters referred to above are taken from the AMBER GAFF force field [67]. Subsequent 2000 step minimizations are performed (a) to relax the water with protein fixed by positional constraints, (b) to relax the protein with all waters constrained, and (c) relax the whole system. Following this minimization protocol, all systems are equilibrated at constant pressure and temperature (NPT) conditions for 1 ns, with the protein fixed by positional constraints. The pressure is regulated using the Berendsen barostat [75] with isotropic position scaling (ntp = 1) and a pressure relaxation time of 2.0 ps, while a Langevin thermostat [76] with collision frequency of 2.0 ps 21 is used to increase the temperature of the system from 0 to 300K. The protein constraints are then lifted and a subsequent 2 ns NPT equilibration is performed at 300K to verify the density of the system is reasonable and stable. The last equilibration step is performed at constant volume and temperature (NVT) for 5 ns at 300K to prepare the system for production MD simulations. All dynamics are conducted using the pmemd.cuda engine [66,77], with Particle Mesh Ewald summations used for computing long-range electrostatic interactions and short-range nonbonded interactions truncated beyond a cutoff of 10 Å [78,79].

Accelerated molecular dynamics (aMD) simulations of IspH
Given current computational power, most MD simulations are limited to sampling timescales on the order of 10-1000 ns. Since many biomolecular processes, including, for example, protein folding, ligand binding, and cis/trans isomerization events, may occur on the order of milliseconds to days, enhanced sampling techniques that facilitate traversing of configuration space efficiently are often implemented to provide information about the relevant conformations of biomolecules [80,81]. Accelerated molecular dynamics (aMD) simulations promote enhanced sampling of systems without the need for defining a reaction coordinate. In aMD simulations, when the potential energy of the system, V(r), is below a threshold energy level, E, a boost energy, DV(r), is applied to encourage exploration of other areas of phase space (Eq. 1). The parameter a modulates the aggressiveness of this boost by altering the depth of the wells in the modified potential.
We employ the dual-boost implementation of aMD to boost both dihedral and total potential energy force field terms to promote side chain dihedral angle rotations and diffusive transitions, respectively [82,83]. We set the parameters E and a for our systems by defining these variables for the dihedral and total potential energy components with respect to the number of residues in the system, N res , and the number of atoms in the system, N atoms , respectively (Eq. 2-5): Subsequent reweighting of the trajectory frames from the aMD simulations using a tenth-order Maclaurin series expansion allows us to extract canonical ensemble averages of the system (further details included in Text S2). Recently, both these methodologies for obtaining aMD parameters and reweighting aMD results were successfully applied to bovine pancreatic trypsin inhibitor (BPTI) to properly obtain the relative populations of relevant, low-lying energetic states [84]. For some semblance of statistics, 36100 ns aMD simulations are performed on all systems explored in this study.

Molecular dynamics analysis
RMSD, RMSF, clustering, and dihedral angle analyses are all performed using the AmberTools 12 suite of programs [66]. Alignment procedures implemented in the RMSD and tRMSF calculations are performed with respect to the [Fe 3 S 4 ] + (open, substrate-free) IspH crystal structure (PDB ID: 3DNF [23]), aligning to the backbone atoms of D1, as this domain is the most rigid in all simulations. Clustering analyses for each of the simulations use pairwise RMSD computed for C a atoms between frames to divide the cumulative trajectories for each system simulated into eight clusters using the average-linkage algorithm [85].  [42,44] processes the dihedral angle distributions for each of these 50 ns blocks as computed by the g_torsion program from the GROMACS suite of programs [91], computing the KL Divergence for a specific dihedral angle using Eq. 6:

Principal component analysis
In this equation, p i refers to the probability that a particular dihedral angle from the [Fe 4 S 4 ] 2+ /HMBPP (closed) simulations falls into a specific range of torsional space, which has been divided into 12u bins. The term p i * is the corresponding probability that the same dihedral angle from the [Fe 4

S 4 ] 2+
(open,substrate-free) simulation falls into the same bin. Combining the KL terms for each of the dihedral angles (Q, y, and x's) of a given residue provides a value for the KL divergence of a specific residue:

Sequence conservation
Using the Bio3D suite of programs, we compute the Shannon entropy [45,46] according to equation 3 for all residues in the A. aeolicus IspH structure with a 22-letter alphabet, where the 20 amino acids are included, as well as a gap character '-' and a mask character 'X' [40]. After normalizing the Shannon entropy score, residues that are fully conserved assume the value 1, while residues with no conservation have a Shannon entropy of 0.

FTMAP
We employ FTMAP [48] to allow many drug-like, organic fragments to bind to representative structures from the dominant clusters from [Fe 4

S 4 ] 2+
(open,substrate-free) aMD simulations. FTMAP utilizes a fast Fourier transform (FFT) algorithm to allow the organic probes to sample many positions along the protein surface. Prevalence of fragment hits along the protein surface signifies ''hot spots'' that correspond to potentially druggable pockets [48].
We measure the ability of residues in substrate-free IspH to bind FTMAP probes by first defining binding of the residue by the probe as existent if the distance between their respective heavy atoms is less than 5 Å . We then combine these binding results across all dominant clusters from the [Fe 4

S 4 ] 2+
(open,substrate-free) aMD simulations and count the number of probes that bind each residue. Normalization of these data indicates the relative propensity of each residue to bind drug-like molecules [92].  Figure S5 Visual representation of the [Fe 4 S 4 (SCH 3 ) 3 OH 2 ] 12 model cluster utilized to obtain charges for the [4Fe-4S] 2+ cluster and its coordinating thiolate residues. Atom labels correspond to those accompanying charges in Table S1. (TIF) Figure S6 Atom labels that correspond to the charges and atom types for the HMBPP molecule given in Table S3. (TIF)   2+ cluster in simulations of IspH (references 12 and 13 in Text S1). (PDF)

Table S3
Force field parameters used for HMBPP. The atom types listed are assigned their respective nonbonded parameters in the AMBER GAFF force field [67]. (PDF) Text S1 Procedure for obtaining force field parameters for the [4Fe-4S] cluster, its coordinating cysteines and HMBPP. (PDF) Text S2 Procedure employed for reweighting aMD trajectories in the principal component analyses. (PDF)