Covalent docking predicts substrates for haloalkanoate dehalogenase superfamily phosphatases.

Enzyme function prediction remains an important open problem. Though structure-based modeling, such as metabolite docking, can identify substrates of some enzymes, it is ill-suited to reactions that progress through a covalent intermediate. Here we investigated the ability of covalent docking to identify substrates that pass through such a covalent intermediate, focusing particularly on the haloalkanoate dehalogenase superfamily. In retrospective assessments, covalent docking recapitulated substrate binding modes of known cocrystal structures and identified experimental substrates from a set of putative phosphorylated metabolites. In comparison, noncovalent docking of high-energy intermediates yielded nonproductive poses. In prospective predictions against seven enzymes, a substrate was identified for five. For one of those cases, a covalent docking prediction, confirmed by empirical screening, and combined with genomic context analysis, suggested the identity of the enzyme that catalyzes the orphan phosphatase reaction in the riboflavin biosynthetic pathway of Bacteroides.

With the explosion of protein sequences, protein functional assignment has emerged as a key problem of the postgenomic era.[1] Despite much progress,[2] sequence-based bioinformatics approaches are mostly limited to annotation transfer of known functions.[3] Meanwhile, function prediction using structure alone is also challenging,[4−6] in part due to the multiple chemical reactions catalyzed by enzymes sharing the same folds.[7−9] Structure-based methods have had most success when they have been combined with ligand chemistry, often via molecular docking.[10−17] In these calculations, libraries of candidate substrates are fit into active sites. Noncovalent complementarity between the protein and the ligand is calculated, using either high-energy intermediate[10,11] or ground-state[18,19] forms of the candidate substrates. Whereas this method suffers from the well-known weaknesses of docking,[20,21] it has nevertheless succeeded in predicting the activities of several families of enzymes, and a much larger number of individual enzymes by annotation transfer.

A key gap in this docking approach has been the reliance on modeling noncovalent fit between a substrate and an enzyme, using modifications of methods first developed for inhibitor discovery.[22−26] Whereas this has proven effective for metalloenzymes such as those in the amidohydrolase and enolase superfamilies, many enzymes proceed through a covalent intermediate that does not lend itself readily to noncovalent modeling. For instance, serine proteases[27] and esterases[28] proceed through an acyl-enzyme intermediate, as do β-lactamases,[29] while decarboxylases and transaminases often form covalent adducts with PLP cofactors.[30] Indeed, some have speculated that many enzymes undergo covalent reactions in the key recognition step along the reaction coordinate.[31] For these enzymes, noncovalent docking of candidate substrates is problematic, as the bond-length approach of the covalent intermediate, and the constraints of the new covalent bond, are poorly modeled by the noncovalent terms of standard docking.

We were thus inspired to investigate the application of a new covalent docking screening method, DOCKovalent,[32] to substrate prediction for enzymes that proceed through covalent intermediates. The method combines covalent bond-length and angle constraints with noncovalent complementarity, drawn from standard docking, and enables large-scale library screens. As with classical, noncovalent docking, the method makes important approximations and adds new ones. Most importantly, it does not calculate the energy of the covalent terms (bond length and angle terms are ignored, as are new torsional energies) but relies exclusively on restraints to model the covalent adduct and complementarity energies from the noncovalent terms. Whereas this has advantages—preventing, for instance, the dominance of covalent terms—the approximation is substantial; as is true with any docking method, it must be tested experimentally before it can be shown to be useful. While covalent docking was used in the past retrospectively to predict substrates of glutathione transferases[33] and predict the chain length of polyprenyl transferases substrates,[34] to our knowledge it was never used in large scale against an enzyme family with a diverse substrate range.

Here we describe the testing of this covalent screening approach against enzymes of the haloalkanoate dehalogenase (HAD) superfamily (HADSF), a superfamily with almost 80 000 sequences in the Structure–Function Linkage Database.[35] Largely dominated by phosphatases, HAD enzymes have wide substrate diversity,[36] with substrates ranging from phosphoserine and histidinol-phosphate through sugar and nucleotide phosphates. We undertake parallel docking and screening campaigns of exactly the same library of phosphate-bearing candidate substrates, using a recently developed empirical screening method to experimentally characterize the substrate specificity of multiple HADSF enzymes. Evaluation of the covalent docking approach by both retrospective enrichment of known substrates and prospective docking for new substrates is considered. Guided by the docking, empirical screens, and by genomic context, we suggest a probable enzyme to fulfill an orphan reaction in the riboflavin biosynthetic pathway, the catalysis of phosphate hydrolysis from 5-amino-6-(5-phospho-d-ribitylamino) uracil. This reaction was hypothesized decades ago, and only recently it was suggested that HADSF members might fill this role[37] in flavinogenic microorganisms.

Methods

Library Generation

The 167 phosphate substrates (Supplementary Dataset 1, Supporting Information) were represented as isomeric smiles strings. The oxygen atom connecting the substrates to the target phosphate was converted to a “dummy” atom to represent the location of the covalent bond. Ligand conformations were generated using Omega[38] as described in ref (32). Corina[39] (Molecular Networks, Erlangen, Germany) was used to generate initial 3D structures and stereoisomers. EPIK[40] (Schrodinger software, Catsville, NY) was used for protonation and tautomer assignment. AMSOL[41] was used to assign partial charges and solvation energies.

Docking

DOCKovalent[32] is a covalent adaptation of DOCK3.6.[42,43] Given a pregenerated set of ligand conformation and a covalent attachment point, it exhaustively samples ligand conformations around the covalent bond and selects the lowest energy pose using a physics-based energy function. For the docking reported in this work, a pentavalent triagonal bipyrimidal phosphate was first modeled covalently attached to the catalytic aspartate residue, either by docking or by manual placement. The substrate library was then covalently docked to the remaining axial phosphate oxygen with a bond length of 1.4 ± 0.3 Å sampled in 0.1 Å increments, bond angle (P–O–ligand) of 120 ± 10° in 2.5° increments, and bond angle (O-ligand covalent attachment point-rest of ligand) of 109.5 ± 10°, also in 2.5° increments. Scoring was as described[32] using a physics-based energy function which uses precalculated van der Waals, electrostatics (calculated with DELPHI[44]), and solvent-excluded desolvation[42] grids, with dampening of the electrostatic potential of the phosphate oxygens to avoid excessive interactions of the ligand with its own phosphate. The receptor is kept fixed throughout the docking simulation.

Retrospective Assessment

For the pose recapitulation benchmark, the phosphate oxygen was used as an anchor point. RMSD was calculated for all substrate heavy atoms excluding the phosphate. The screening hit-rate was defined as the number of substrates divided by the size of the tested library. The docking hit-rate was defined as the number of substrates ranking in the top 20 of the docking hit list divided by 20. Docking enrichment was calculated by dividing the docking hit rate by the screening hit rate; random selection gives an enrichment of 1.

Noncovalent Docking of High-Energy Intermediates

Here the phosphate reactive center is also represented in a trigonal bipyramidal geometry. Three-dimensional models of the ligands are prepared with Corina,[39] and the molecules are then transformed to their intermediate state using OEChem TK, version 1.7.4, OpenEye Scientific Software, Inc., Santa Fe, NM, USA, www.eyesopen.com, 2010. In these reactions, first the phosphate is attacked with a hydroxide and second the phosphate protonation states are enumerated, creating multiple high-energy intermediates for each molecule. Next the molecules are prepared as described previously.[10,45]

Receptor structures were prepared for noncovalent docking as described.[10,45,46] Chain A from each structure was chosen, except for 3DDH where chain B was used. 1RKU, 3PGV, 3N07, 2OBB, 3GYG, 3MMZ, and 3N1U were prepared as dimers. Only the highest occupancy rotamer of each residue was kept, and incomplete side chains were replaced by using rotamer libraries.[47] The charge of the Mg2+ ion was reduced to +1.4 and two water molecules were placed to coordinate the metal ion in the same idealized geometry for all structures. Parallel docking experiments were run with the wild type enzyme, and with a mutation of the catalytic aspartate to a serine residue, in which the hydroxyl was oriented to coordinate the Mg2+. The high-energy library was then docked to each active site with DOCK3.6.[42] The docking was run using receptor and ligand bin sizes of 0.4 Å, an overlap of 0.1 Å, a distance tolerance of 1.5 Å, color matching turned off, and 250 cycles of rigid-body minimization.

Genome Neighborhood Networks

Beginning with gi|29346380, a set of similar protein sequences (seed sequences) were collected by performing a BLAST[48] search against the NCBI GenPept database[49] at an e-value cutoff of 10–20. For each protein BLAST hit, the analogous gene was determined, and the proteins corresponding to the 10 genes upstream and downstream to each were collected. All of these proteins were combined, and an all-by-all BLAST was performed at an E-value cutoff of 10–14 against a custom database containing only the proteins of interest. A cytoscape[50] network was created from these BLAST results. Each node in the network represents a single protein sequence, and each edge represents the pairwise connection with the most significant BLAST E-value (better than the cutoff) connecting the two sequences. Connections between nodes are only shown if the E-value of the best Blast hit between two sequences is at least as good as the specified E-value cutoff. Lengths of edges are not meaningful except that sequences in tightly clustered groups are more similar to each other than sequences with few connections. The nodes were arranged using the yFiles organic layout provided with Cytoscape version 2.8. Annotation information retrieved from Swiss-Prot[51] (functional annotation) and NCBI[49] (taxonomy information; position of the corresponding gene relative to the seed gene) was associated with each node as applicable.

Empirical Screen

The target protein was diluted to 10 μM in a volume of 2.5 μL in assay buffer (20 mM HEPES, pH 7.5, 100 mM NaCl, 5 mM MgCl2, 2 mM DTT) and added to each of 167 putative substrates diluted in assay buffer to 2 mM, in a volume of 2.5 μL (see Supplementary Dataset 1, Supporting Information for the substrate list) in duplicate wells using 384-well plates (Corning 384 Well low Volume Black Clear Bottom - Cat. No. 3540) and incubated for 30 min followed by addition of BioMol Green (13.5 μL; Prod. No. BML-AK111), a dye that is sensitive to the presence of free phosphate. The mixture was incubated for an additional 45 min and absorbance at 650 nm was read using an EnVision multilabel reader (product number: 2104-0010).

Preparation of Recombinant Escherichia coli RibA

The DNA encoding the ribA gene from E. coli was amplified by polymerase chain reaction (PCR) using E. coli genomic DNA (ATCC 29148D), Pfu Turbo DNA polymerase, and oligonucleotide primers (5′-AATCATATGCAGCTTAAACGTGTG and 5′-GGTTATTTTGGATCCGGCAAGC) containing restriction endonuclease cleavage sites NdeI and BamHI. The pET-15b TEV vector, cut by restriction enzymes NdeI and BamHI, was ligated to the PCR product that had been purified and digested with the same restriction enzymes. The ligation product was used to transform to NEB Express Iq Competent E. coli cells that were then grown on a ampicillin-containing agar plate. A selected colony was checked for RibA expression, and the isolated plasmid was sequenced to verify the correct gene sequence. For RibA preparation, the transformed cells were grown at 37 °C with agitation at 225 rpm in 1 L Terrific Broth containing 100 μg/mL ampicillin to an OD600 of 0.6 and then induced for 12 h at 16 °C with 0.4 mM isopropyl beta-d-thiogalactopyranoside (IPTG). The cells were harvested by centrifugation (6500 rpm for 15 min at 4 °C) to yield 16 g/L of culture medium. The cell pellet was suspended (1 g of wet cells/10 mL) in ice-cold buffer A (20 mM HEPES (pH 7.5), 5 mM MgCl2, and 10 mM imidazole). The cell suspension was passed through a French press at 1200 psi before centrifugation at 20 000 rpm at 4 °C for 20 min. The supernatant was loaded onto 5 mL Ni-NTA agarose column at 4 °C. After the column had been washed with 100 mL of buffer B (20 mM HEPES (pH 7.5), 300 mM NaCl, and 20 mM imidazole), the enzyme was eluted with 200 mL of elution buffer C [20 mM HEPES (pH 7.5), 50 mM NaCl, and 500 mM imidazole]. The column fractions were analyzed by SDS–PAGE, and the desired fractions were combined before dialysis at 4 °C against 6 L OF triethanolamine hydrochloride (TEA) buffer (20 mM TEA (7.5) and 100 mM NaCl). The final yield was 20 mg of RibA/g of wet cells.

Preparation of Recombinant Thermotoga maritima RibD

The plasmid (PSI biology) containing the T. maritima ribD gene was transformed to BL21(DE3)pLysS competent cells and then grown on an ampicillin-containing agar plate. A selected colony was checked for RibD expression. For RibD preparation, the transformed cells were grown at 37 °C with agitation at 225 rpm in 2 L of Terrific Broth containing 100 μg/mL ampicillin to an OD600 of 1.0 and then induced for 12 h at 25 °C with 0.6 mM IPTG. The cells were harvested by centrifugation (6500 rpm for 15 min at 4 °C) to yield 12 g/L of culture medium. The cell pellet was suspended (1 g of wet cells/10 mL) in ice-cold buffer A (50 mM HEPES (pH 8.0), 50 mM NaCl, and 10 mM imidazole). The cell suspension was passed through a French press at 1200 psi before centrifugation at 20 000 rpm at 4 °C for 20 min. The supernatant was loaded onto 5 mL Ni-NTA agarose column at 4 °C. After the column had been washed with 100 mL of buffer B (50 mM HEPES (pH 8.0), 300 mM NaCl, and 40 mM imidazole, 10% glycerol), the enzyme was eluted with 200 mL of elution buffer C (50 mM HEPES (pH 8.0), 300 mM imidazole, 10% glycerol). The column fractions were analyzed by SDS–PAGE, and the desired fractions were combined and then dialyzed at 4 °C against 6 L 50 mM HEPES (pH 8.0), 50 mM NaCl. The final yield was 4.9 mg of RibD/g of wet cells.

Preparation of Recombinant E. coli Putative 5-Amino-6-(5-phospho-d-ribitylamino) Uracil Phosphatase

The cell stock [BL21 (Ros2), EFI:501083, UniProt: Q8A947] was grown at 37 °C with agitation at 225 rpm in 10 mL of Terrific Broth containing 100 μg/mL ampicillin overnight, and then the 10 mL culture was transferred to 2 L of Terrific Broth containing 100 μg/mL ampicillin and grown at 37 °C with agitation at 225 rpm to an OD600 of 0.8, and then induced for 12 h at 16 °C with 0.5 mM IPTG. The cells were harvested by centrifugation (6500 rpm for 15 min at 4 °C) to yield 9 g/L of culture medium. The cell pellet was suspended (1 g of wet cells/10 mL) in ice-cold buffer A (25 mM HEPES (pH 7.5), 50 mM NaCl, and 5 mM MgCl2, 5% glycerol, 10 mM imidazole). The cell suspension was passed through a French press at 1200 psi before centrifugation at 20 000 rpm at 4 °C for 20 min. The supernatant was loaded onto 5 mL Ni-NTA agarose column at 4 °C. After the column had been washed with 100 mL of buffer B (25 mM HEPES (pH 7.5), 300 mM NaCl, 5 mM MgCl2 and 40 mM imidazole), the enzyme was eluted with 200 mL of buffer C (25 mM HEPES (pH 7.5), 5 mM MgCl2, 500 mM imidazole). The column fractions were analyzed by SDS–PAGE, and the desired fractions were combined and then dialyzed at 4 °C against 6 L 25 mM HEPES (pH 7.5), 50 mM NaCl, 5 mM MgCl2, 5% glycerol. The final yield was 8.6 mg of protein/g of wet cells.

Enzymatic Synthesis of 5-Amino-6-(5-phospho-d-ribitylamino) Uracil

A reaction mixture containing 50 mM HEPES (pH 7.9), 5 mM MgCl2, 4.05 mM GTP, 4.54 mM NADPH, 90 μM RibA, and 175 μM RibD in a final volume of 222 μL was incubated at 37 °C for 1.5 h. The concentration of 5-amino-6-(5-phospho-d-ribitylamino) uracil was determined by the decrease in absorbance at 340 nm, indicative of NADPH oxidation (ε340 = 6220 M–1 cm–1).

Determination of Steady-State Kinetic Constants

Initial velocities for putative 5-amino-6-(5-phospho-d-ribitylamino) uracil phosphatase-catalyzed hydrolysis of sugar phosphate were measured at 25 °C using assay solutions that contained 1 mM MgCl2, 0.1 mM sodium azide 1.0 unit/mL purine nucleoside phosphorylase, and 0.2 mM MESG in 50 mM Tris-HCl (pH 7.5). Absorbance changes were monitored at 360 nm (Δε = 9.8 mM–1 cm–1). The steady-state kinetic parameters (Km and kcat) were determined by fitting the initial velocity data measured at varying substrate concentrations (ranging from 0.5 Km to 5 Km) to the equation: V0 = (Vmax[S])/([S] + Km) where V0 is the initial velocity, Vmax the maximum velocity, [S] the substrate concentration, and Km the Michaelis constant for the substrate, using the SigmaPlot Enzyme Kinetics Module. The kcat values were calculated from Vmax and [E] according to the equation kcat = Vmax/[E], where [E] is the enzyme concentration.

Results

Overview of the Method

Catalysis by HAD phosphohydrolases proceeds via an aspartylphosphate intermediate[52] (Supplementary Figure 1, Supporting Information). Since the phosphate-binding mode is conserved in the HAD catalytic core domain[36] we begin by covalently attaching a pentavalent, trigonal-bipyrimidal phosphate[53] to the catalytic aspartate (Oδ–P bond length = 1.9 ± 0.3Å; Cγ–Oδ–P angle = 120 ± 10°; Figure 1). This can be done manually or by applying the covalent docking protocol described below to a model phosphorylated molecule. Once the pentavalent phosphate is placed, we covalently dock a library of candidate substrates to the axial phosphate oxygen (Figure 1; See Methods for geometric parameters). The substrate library of 167 phosphorylated molecules mirrors that of the library used in the empirical screens (Supplementary Data set 1, Supporting Information). For each ligand, we sample ligand conformations with respect to the covalent bond to the fixed phosphate, constrained by ideal bond lengths and angles, in 20° increments around the two newly formed torsion angles, with conformations for the remaining portions of the molecule precalculated. Each sampled conformation is scored using a physics-based scoring function in DOCK3.6,[42] which evaluates the ligand van der Waals and electrostatic interactions, and corrects for ligand desolvation. The best scoring pose for each ligand is saved. The library is then ranked by this docking score.

Figure 1

Overview of covalent docking to the HADSF. An illustration of the sampling in the HADSF covalent docking. A pentavalent trigonal-bipyramidal phosphate is modeled as covalently attached adduct to the catalytic aspartate, either manually or via covalent docking of the phosphate. Substrates are then covalently docked to the phosphate axial oxygen, exhaustively sampling the two indicated dihedral angles as well as pregenerated ligand conformations. The best pose is kept for each substrate and substrates are than ranked based on their scores.

Inadequacy of Noncovalent Docking

Previously, we used noncovalent docking of high-energy intermediates to predict substrate recognition by amidohydrolase family enzymes,[10−13] which do not proceed via a covalent intermediate. To investigate whether this approach could model enzymes like those of the HAD superfamily, which do proceed by such an intermediate, we docked high-energy intermediates of the phosphatase reaction against both the unmodified enzyme structure and against an artificial form of that structure truncated at the nucleophilic aspartate (substituting the aspartate with a serine) so as to allow close approach of the substrate and enzyme.

While in several of the cases substrates indeed ranked in the top of the docking hit list (Supplementary Table 1, Supporting Information), it was often the case that the predicted docking pose was not competent for the reaction (Supplementary Figure 2, Supporting Information). This supported efforts at modeling the reaction using the covalent docking approach, which by design produces only competent poses.

Retrospective Assessment of Covalent Docking

We first investigated the ability of the covalent docking method to recapitulate the crystallographic poses of HAD–substrate complexes. We assembled a small benchmark of seven liganded structures covering a diverse range of substrates and HAD subtypes (Supplementary Table 2). In five of the seven cases covalent docking recapitulated the substrate binding pose to less than 2Å RMSD (Figure 2), while in the other two structures the overall binding pose was recapitulated but not with atomic accuracy (RMSD = 2.14 Å and 3.1 Å). A more stringent test is docking to unliganded structures. When using ligand free structures (or in one case, a complex with an alternative ligand bound) performance decreased, although still in three cases poses were recovered to less than 2 Å RMSD (Supplementary Table 2).

Figure 2

Substrate pose recovery by covalent docking. Examples of covalent docking pose predictions (magenta) for known HAD (white)/substrate (yellow) complexes. (a) Deoxyribonucleotidase in complex with deoxyuridine (PDB: 2I7D); (b) sucrose-phosphatase in complex with sucrose-6P; (c) human pyridoxal phosphate phosphatase in complex with pyridoxal-5P.

To investigate the ability to find HAD substrates using unliganded structures, we calculated docking structures, in covalent high-energy intermediate form, of the library of 167 candidate phosphate substrates that was used in empirical screens for HAD function. We compared docking performance to that of the empirical screen on a set of 20 HADSF members with a solved unliganded crystal structure (Table 1). For each enzyme we considered the top 20 predictions, out of 167 possible, as putative substrates. We considered two alternatives as true substrates from the empirical screen. The poor definition refers to any molecule that showed a signal larger than 0.1 absorbance units in the empirical screen. The good definition refers to substrates that showed a signal larger than one standard deviation from the mean of all poor substrates signal. While poor substrates are turned over by the enzyme, the cellular function of the enzyme is more likely related to good substrates.

Table 1

Retrospective Assessment of HAD Substrate Prediction

	enrichmenta	examples of correctly predicted substrates (docking rank)b
PDB	good	poor	good	poor
2b82c	3.30	1.98	IMP (17)	GMP (10)
				dUMP (16)
4dcc	3.30	1.89	riboflavin-5-phosphate-(FMN) (9)	l-sorbose-1-phosphate (5)
				d-iditol-6-phosphate (11)
				arabinose-5-phosphate (19)
4dfd	3.30	1.89	riboflavin-5-phosphate-(FMN) (4)	l-sorbose-1-phosphate (9)
				d-allitol-3-phosphate(15)
				d-galactitol-6-phosphate(17)
3d6j	3.30	2.50	d-tagatose-6-phosphate (4)	arabinose-5-phosphate (1)
			2-deoxyribose-5-phosphate (5)	d-ribose-5-phosphate(10)
			d-2-deoxy-ribitol-5-phosphate (6)	erythrose-4-phosphate (13)
1nrw	2.93	2.16	isoerythritol-4-phosphate (7)	O-phosphorylethanolamine (3)
			2-deoxyribose-5-phosphate (18)	d-ribulose-5-phosphate (4)
			glycerol-3-phosphate (19)	d-threitol-4-phosphate (5)
			ribitol-5-phosphate (20)	glycerol-phosphate-(GP) (6)
1rku	2.40	0.60	O-phospho-l-serine (7)	arabinose-5-phosphate (16)
1te2	2.10	2.15	mannose-6-phosphate (11)	glucosamine-6-phosphate (2)
			ribitol-5-phosphate (14)	d-2-keto-glucose-6-phosphate (5)
			2-deoxy-d-glucose-6-phosphate (18)	l-ribitol-5-P (6)
3pgv(AC)	1.55	1.52	glucosamine-6-phosphate (2)	allose-6-phosphate (8)
			d-sedoheptulose-7-phosphate (6)	IMP (9)
			d-psicose-6-phosphate (13)	dGMP (12)
			CMP (20)	d-allitol-6-phosphate (14)
4eek	1.20	1.13	UMP (6)	d-altronate-6-phosphate (14)
			d-allonate-6-phosphate (12)	CMP (16)
				d-allitol-6-phosphate (17)
3r4c	1.20	0.83	dAMP (20)	CMP (7)
				UMP (10)
				arabinose-5-phosphate (17)
3n07(B)	1.20	0.27	d-fructose-1,6-bisphosphate (17)
3niw	0.80	0.44	arabinose-5-phosphate (7)	mannose-6-phosphate (11)
2obb	0.00	2.20		dl-glyceraldehyde-3-phosphate (16)
2hx1	0.00	1.47		d-3-deoxy-glucose-6-phosphate (14)
				mannitol-6-phosphate (17)
3gygd	0.00	1.47		dl-glyceraldehyde-3-phosphate (3)
				d-ribulose-5-phosphate (7)
				l-sorbose-1-phosphate (9)
3s6j	0.00	1.44		erythrose-4-phosphate (2)
				dl-glyceraldehyde-3-phosphate (3)
				arabinose-5-phosphate (6)
3mmzd	0.00	0.69		d-fructose-1,6-bisphosphate (2)
				6-phosphogluconic-acid (5)
				mannose-6-phosphate (14)
3n1ud	0.00	0.65		6-phosphogluconic-acid (4)
				α-d-glucose-1,6-bisphosphate (13)
				arabinose-5-phosphate (14)
1z5g	0.00	0.00
3ddh(B)	0.00	0.00

Docking enrichment calculated separately for good and poor substrates as defined in the text. An enrichment value of 1 corresponds to random prediction. Values larger than 1 indicate successful docking predictions.

Examples of predicted substrates (docking rank indicated in parentheses) that were empirically shown to serve as good or poor substrates for the indicated enzyme.

Unless indicated otherwise chain A was used for prediction.

Crystallographic symmetry was applied to model the “biological” binding site.

In 12 out of the 20 cases, we find one or more good substrate among the top 20 docking predictions, in 11 of these we also observed enrichment of these good substrates 20–230% better than expected by random (Table 1). In six other cases, a poor substrate is detected by docking (four of them with substantial enrichment), and in two cases no substrate is detected by docking. Although these enrichments are substantially better than a random predictor, the enrichment reached statistical significance (p-value < 0.05, as evaluated by Fisher’s exact test; Supplementary Table 3) in only four of these cases, possibly due to the large number of substrates for many of the HAD phosphatases.

Prospective Substrate Prediction

Encouraged by the covalent docking performance in recovering the geometries of bound substrates, and enriching for known vs decoy substrates in retrospective screens, we next covalently docked prospectively against structures of HADSF members for which no function was known. We docked the library of 167 phosphorylated metabolites against the structures of seven such enzymes (Table 2) and in parallel screened the same library against them empirically. Covalent docking was able to predict good substrates in five cases and poor substrates in six cases, with a better than random enrichments in 4/7 cases (Table 2; full empirical screening results are in Supplementary Dataset 1). For instance, d-tagatose-6P ranked 11 by docking for the orphan HAD enzyme EFI-508415 (PDB: 4gxt), forming extensive polar contacts with Asp52, Arg339, and Asp 317 (Supplementary Figure 3a). It was shown by the empirical screen to be a good substrate for this enzyme. l-Xylose-5-P ranked second by docking for the orphan HAD enzyme EFI-502344 (PDB: 3dv9) based on good binding site complementarity and hydrogen bonding network with His38, Trp42, Asp28 and the backbone of Gly129 (Supplementary Figure 3b). Empirically it was observed to be a robust substrate. d-Ribitol-5P ranked 10 by docking to the orphan EFI-501083 (PDB: 4jb3, UniProt: Q8A947), displaying good geometric fit supplemented with two hydrogen bonds to Glu62 (Supplementary Figure 3c). In the empirical screen, it was the best substrate for this enzyme. We decided to further investigate this enzyme, for which the docking prediction combined with genomic context analysis hinted at an interesting biological function.

Table 2

Prospective Prediction of HAD Substrates via Covalent Docking

	enrichmenta	examples of correctly predicted substrates (docking rank)b
PDB	good	poor	good	poor
1nf2c	2.09	1.19	2-deoxy-d-glucose-6-phosphate (14)	acetyl-phosphate (5)
				2-deoxyribose-5-phosphate (13)
4gxt	1.79	1.28	d-tagatose-6-phosphate (11)	6-phosphogluconic-acid (2)
			ribitol-5-phosphate (12)	l-sorbose-1-phosphate (6)
			mannitol-6-phosphate (20)	d-glycero-beta-d-manno-heptose-1,7-bisphosphate (7)
4jb3	1.67	1.34	ribitol-5-phosphate (10)	l-ribose-5-P (6)
			meso-erythritol-4-phosphate (14)	l-ribitol-5-P (11)
				l-lyxitol-5-P (13)
3dv9	0.70	1.92	d-threitol-4-phosphate (18)	d-xylose-5-P (1)
				l-xylose-5-P (2)
				d-2-keto-glucose-6-phosphate (3)
2b0c	0.38	0.97	α-d-glucose-1-phosphate (16)	allose-6-phosphate (3)
				d-tagatose-6-phosphate (5)
				d-glucose-6-phosphate (6)
2fi1	0.00	0.60		dTTP (13)
1ydf	0.00	0.00

Docking enrichment calculated separately for good and poor substrates as defined in the text. An enrichment value of 1 corresponds to random prediction. Values larger than 1 indicate successful docking predictions.

Examples of predicted substrates (docking rank indicated in parentheses) that were empirically shown to serve as good or poor substrates for the indicated enzyme. See Supplementary Data set 1 for additional substrates discovered for these targets.

Chain A was used for prediction for all targets.

EFI-501083 Catalyzes an Orphan Reaction in the Riboflavin Biosynthesis Pathway

EFI-501083 (PDB: 4jb3, UniProt: Q8A947) is a HAD enzyme from Bacteroides thetaiotaomicron with unknown function, whose structure was determined as part of a larger effort to determine enzyme function.[54] Empirical screens against this enzyme revealed it was a broad range 5–6 carbon alcohol and aldol sugar phosphatase (Supplementary Data set 1). Consistent with this observation, 5 of the top 20 prospective docking predictions were confirmed as substrates by the empirical screen, with the best substrate, as mentioned, d-ribitol-5-phosphate, ranked 10 by docking (Supplementary Figure 3c).

Bacterial enzymatic pathways are often organized in operons such that enzymes of the same pathway are found in sequence proximity in the genome. For this reason the genomic context, here defined as 10 genes upstream and downstream of the target gene, can provide important clues for an enzyme’s function.[15] Genome neighborhood networks[55] expand this notion and incorporate the genomic context of closely related sequences (see Methods).

The enzyme EFI-501083 (UniProt: Q8A947) was used to seed a genome neighborhood network calculation. One of the closest clusters in the average genomic neighborhood contained enzymes annotated as participating in the riboflavin biosynthesis pathway (Supplementary Figure 4). Specifically, one cluster contained RibF, a bifunctional riboflavin kinase, and FMN adenylyltransferase, which is comprised of two enzymes (fused in most bacteria) of the common FMN/FAD pathway.[56,57] Another cluster contained RibD, a reductase that catalyzes the conversion of 5-amino-6-(5-phospho-d-ribosylamino) uracil (compound 1) to 5-amino-6-(5-phospho-d-ribitylamino) uracil (compound 2; Figure 3a). The next step in riboflavin biosynthesis is the removal the phosphoryl group from compound 2 by a putative, unidentified, phosphatase.[58,59] This genomic context hint, combined with the docking prediction and screening evidence that 5-phospho-d-ribitol, a substructure of 2 (Figure 3a; marked in red), is a substrate for EFI-501083, led us to suspect that it is in fact the missing phosphatase of this pathway in B. thetaiotaomicron.

Figure 3

EFI-501083 (UniProt: Q8A947) may be the missing enzyme for the orphan reaction in the riboflavin biosynthetic pathway. (a) Two consecutive intermediate steps in the riboflavin biosynthetic pathway.[58] The first is a reduction catalyzed by RibD, which is found in the genome neighborhood network of EFI-501083. The second is a phosphohydrolase reaction that may be catalyzed by EFI-501083. Ribitol-5P — a substructure of 2 (marked in red) was predicted by docking as a substrate for this enzyme and was a good substrate discovered by empirical screening. (b) A model of 2 in the binding pocket of EFI-501083 suggests it can indeed bind the full substrate. (c) Phosphatase catalytic activity of EFI-501083 on substrate 2.

While docking 2 to the structure (PDB: 4jb3) did not provide a fit in the enzyme binding pocket, repacking of the side-chains of Phe83 and Ile16 and subsequent minimization resulted in an enzyme structure that could accommodate 2 in the docking pose of d-ribitol-5-phosphate (Figure 3b). This prompted us to assess 2 as a substrate for EFI-501083 in vitro.

In a phosphatase assay against EFI-501083, d-ribitol-5-phosphate and the related d-ribitol-1-phosphate were found to be enzyme substrates, with kcat/KM of 1.4 × 103 M–1 s–1 and 1.7 × 103 M–1 s–1, respectively (Table 3). Next we measured the kinetics of phosphoryl hydrolysis of 2 generated by enzymatic synthesis using RibD.[60] The resulting kcat/KM of 1.6 × 102 M–1 s–1 (Table 3; Figure 3c) is within 10-fold of the rates measured for the preceding step in the riboflavin pathway[60] consistent with the hypothesis that EFI-501083 may serve as the missing phosphatase of this biosynthesis pathway.

Table 3

Steady-State Kinetic Constants for Putative 5-Amino-6-(5-phospho-d-ribitylamino) Uracil Phosphatasea

substrate	kcat (s–1)	Km (μM)	kcat/Km (M–1 s–1)
d-ribitol-1-Pb	1.8 ± 0.1	1030 ± 180	1.7 × 103
d-ribitol-5-Pb	1.5 ± 0.1	1080 ± 140	1.4 × 103
5-amino-6-(5-phospho-d-ribitylamino)uracilc	0.016 ± 0.001	100 ± 10	1.6 × 102

Catalyzed hydrolysis of sugar phosphate at 25 °C and pH 7.5 (see Methods).

Substrate concentration were determined by full conversion of free phosphate catalyzed by EFI-501083.

Substrate concentration was determined by decrease in absorbance at 340 nm, indicative of NADPH oxidation (ε340 = 6220 M–1 cm–1) during the enzymatic synthesis.

Discussion

Two key observations emerge from this study. First, covalent docking captures the substrate recognition for multiple subtypes of HAD superfamily enzymes, which have a key covalent intermediate along their reaction coordinate; this approach should be applicable to enzymes from other superfamilies that go through such an intermediate. Such substrate recognition, both retrospective and prospective, was better captured by the new covalent docking approach than by classical, noncovalent docking. Second, when combined with genome neighborhood networks, the method illuminated the identity of a long hypothesized phosphatase that completes the riboflavin biosynthesis pathway.

Noncovalent docking was proven successful in various modeling and substrate prediction applications[11−13] including those performed on dehalogenases.[61,62] However, it requires no great leap to imagine that classical, noncovalent docking will struggle to model the covalent intermediates that feature in so many enzyme reaction coordinates. For these reactions, the enzyme active site is preorganized to stabilize such a covalent intermediate;[63,64] hence any noncovalent positioning of the intermediate will be suboptimal. Naturally, just because noncovalent docking might be expected to struggle does not mean that a covalent docking method will succeed. The ability to accurately recover substrate-bound crystallographic structures, even using unliganded receptor structures in some cases, and to identify true substrates both retrospectively and prospectively, supports the use of this method, and likely many covalent docking approaches, for prospective substrate prediction. For instance, the accuracy of pose recapitulation for known cocrystal complexes (Figure 2; Supplementary Table 2) leads us to believe that highly scoring covalent docking poses are in productive, on-path conformations, in contrast to the incompatible poses predicted by noncovalent docking (Supplementary Figure 2). The correct prediction of good substrates for 17 out of 27 overall different proteins, five of them prospectively, and poorer substrates in seven additional cases, speaks to the pragmatic applicability of this approach.

The discovery that EFI-501083 can fulfill the role of the long-hypothesized phosphatase in the riboflavin biosynthesis pathway illuminates the potential of this method for new biological discovery. Since the 1950s, when the observation was made that the yield of vitamin B2 could be increased during fermentation of Eremothecium ashbyii by the addition of purines to the culture medium,[65] numerous studies have demonstrated that the atoms of the purine ring system are incorporated into riboflavin and have elucidated the related biochemical transformations (for review see ref (59)). However, although the intermediacy of the pyrimidine 2 is well established, as is the use of the dephosphorylation product 3 by lumazine synthase, the identity of the phosphatase itself has remained elusive. Recently, a clue to the involvement of HADSF members in this phosphatase activity came from plants, which use the same pathway as eubacteria, in which a HADSF member purified from chloroplasts was shown to have FMN hydrolase activity.[66] Newly published work in E. coli presented a chemi-enzymatic synthesis of 2 and evidence that both the HADSF members YigB and YbjI (17% sequence identity) can support turnover of 2, with Km values in the physiological range, 20 and 70 μM, respectively, and specific activity similar to other enzymes in the riboflavin pathway.[37] The single deletion strains of the yigB or ybjI genes grow normally in the absence of riboflavin, a finding that is ascribed to the fact that either can carry out the phosphatase activity and on the observed activity on FMN of other E. coli HADSF enzymes.[67] Here, the candidate substrates prioritized by covalent docking, together with empirical screens and genome context, identified enzyme EFI-501083 from B. thetaiotaomicron as the long sought phosphatase. This enzymatic function could not have been inferred by using each of these prediction methods separately, nor was this function clear from sequence identity alone, as an isofunctional orthologue of either yigB or ybjI, given the sequence identities of only 17% and 5%, respectively. We note that in B. thetaiotaomicron, a second member EFI-501088 with 33% identity to EFI-501083, has activity against FMN[68] with similar low identity to yigB or ybjI (sequence identity 17% and 16%, respectively).

While an advantage of covalent over noncovalent docking is its dramatic reduction in the number of sampled degrees of freedom, it shares the limitations of most docking approaches and must overcome several challenges. We do not pretend here to have met all of these challenges here —key problems such as accurate accounting for the covalent-bond energetics and binding-site flexibility remain unresolved. HAD enzymes can be found in an “open” or “closed” conformation. While for this study we focused on structures solved in the closed conformation, being able to model this conformational change would increase the coverage of putative targets for our approach. Finally, the approach, like docking in general, suffers from false positive predictions. Indeed, compared to its application to inhibitor discovery,[32] which was able to discover potent, reversible, and highly ligand efficient inhibitors of several enzymes with high hit rates, the covalent docking method often only had modest hit rates with many false negatives. We suspect that this reflects the challenges of recognizing phosphorylated ligands more than intrinsic challenges with substrate versus inhibitor prediction, but this remains an area of ongoing research.

These limitations should not obscure the main observations of this study: against 17 different HAD enzymes, covalent docking captured the essential features of substrate recognition, both retrospectively and, more compellingly, prospectively. The prediction that ribitol-5P is a substrate of EFI-501083, combined with genomic context and kinetic analysis implicates it as the enzyme that catalyzes the orphan phosphatase reaction in the riboflavin biosynthetic pathway sought for over 50 years, illustrating the potential applications of the method to pathway discovery and deeper biology.