Identification of Fatty Acid Binding Protein 5 Inhibitors Through Similarity-Based Screening.

Fatty acid binding protein 5 (FABP5) is a promising target for development of inhibitors to help control pain and inflammation. In this work, computer-based docking (DOCK6 program) was employed to screen ∼2 M commercially available compounds to FABP5 based on an X-ray structure complexed with the small molecule inhibitor SBFI-26 previously identified by our group (also through virtual screening). The goal was discovery of additional chemotypes. The screen resulted in the purchase of 78 candidates, which led to the identification of a new inhibitor scaffold (STK-0) with micromolar affinity and apparent selectivity for FABP5 over FABP3. A second similarity-based screen resulted in three additional hits (STK-15, STK-21, STK-22) from which preliminary SAR could be derived. Notably, STK-15 showed comparable activity to the SBFI-26 reference under the same assay conditions (1.40 vs 0.86 μM). Additional molecular dynamics simulations, free energy calculations, and structural analysis (starting from DOCK-generated poses) revealed that R enantiomers (dihydropyrrole scaffold) of STK-15 and STK-22 have a more optimal composition of functional groups to facilitate additional H-bonds with Arg109 of FABP5. This observation suggests enantiomerically pure compounds could show enhanced activity. Overall, our study highlights the utility of using similarity-based screening methods to discover new inhibitor chemotypes, and the identified FABP5 hits provide a strong starting point for future efforts geared to improve activity.

Fatty acid binding proteins (FABP) are a family of lipid chaperone proteins that transport fatty acids and other lipophilic substances. The 10 mammalian FABP isoforms are widely expressed in humans, each with distinct tissue expression patterns and ligand binding preferences.[1,2] Although the amino acid sequence identity of different FABPs range from ∼20% to ∼70%, they share highly similar tertiary structure and binding site conformation.[2] A wide range of physiological functions have been proposed and studied for FABPs. In general, they are involved in transporting fatty acids and other lipophilic ligands to various intracellular sites for metabolism, storage, and signal transduction.[1−3]

Epidermal-type fatty acid binding protein (FABP5, E-FABP), in particular, has been demonstrated to be involved in the N-acylethanolamine (NAE) regulation pathways.[3] NAEs are a family of signaling lipids including the endocannabinoid anandamide (AEA) that activates cannabinoid receptors (CB) and oleoylethanolamide (OEA) and palmitoylethanolamide (PEA) that largely signal through nuclear peroxisome proliferator-activated receptor alpha (PPARα).[4−6] It has been demonstrated through various experiments that visceral, inflammatory, and neuropathic pain can be alleviated by inhibiting fatty acid-amide hydrolase (FAAH), the principal NAE hydrolyzing enzyme, with the subsequent activation of CB and PPARα receptors.[7−10] FABP5 acts as an intracellular shuttle to bring the hydrophobic NAEs through the aqueous environment of the cytoplasm to FAAH for catabolism.[3,11] Pharmacological inhibition or genetic elimination of FABP5 results in significantly elevated NAE levels, thereby increasing activation of CB and PPARα receptors, resulting in antinociceptive and anti-inflammatory effects.[12−14] Consequently, there is evidence that development of small molecule FABP5 inhibitors is a viable path toward a new class of analgesic therapeutics. However, knockout of the closely related heart-type fatty acid binding protein (FABP3, H-FABP) has been implicated in causing cardiac hypertrophy in rodent models,[15] therefore compounds selective against the FABP3 isoform are desirable.

As a potentially important pharmacological target for pain, inflammation, and amelioration of drug withdrawal effects (as seen in rodents), efforts by a number of groups, including our own, have been directed into developing effective small molecule FABP inhibitors.[1,12,13,16−21] Due to the lack of a viable FABP5/small molecule cocomplex at the time, prior work by our group as reported in Berger et al.,[12] employed large-scale virtual screening (DOCK program) to identify FABP leads using an X-ray structure of FABP7 complexed with oleic acid (with oleic acid removed). Importantly, the campaign identified several active compounds represented by SBFI-26 (Ki = 0.93 ± 0.08 μM)[12] with ∼3-fold selectivity for FABP5[20] against FABP3. Follow up work included the design, synthesis, and testing of a large series of SBFI-26 analogues, although none of the analogues showed significant improvements in both affinity and selectivity for FABP5.[13,20] Subsequent X-ray structures of (S)-SBFI-26 complexed with FABP5 and FABP7 provided additional structural insights on small molecule inhibitor binding to FABPs. Notably, the DOCK-predicted binding geometry for SBFI-26 (termed binding pose) with FABP7 corroborated key aspects of the experimentally determined pose despite a difference in ligand stereochemistry.[12,19]

The goal of the current study is to expand upon our previous FABP5 inhibitor development utilizing the newly solved crystal structure of FABP5 cocomplexed with SBFI-26 and a much larger library for virtual screening. The main objectives were 4-fold: (1) perform an in silico screening of ∼2 million drug like small molecules from the ZINC database[22] and prioritize docked poses using similarity-based scoring functions in DOCK using the X-ray pose of SBFI-26 as a reference, (2) purchase and experimentally test a diverse subset of the most promising compounds, (3) prioritize and purchase additional analogues based on any experimentally identified hits for a second iteration of experimental testing, and (4) conduct a comprehensive structural and energetical analysis on compounds showing experimental activity. As detailed below, this approach has led to the identification of a series of nonacid compounds that inhibit FABP5 with moderately strong affinities.

Methods

Structure Selection

The crystallographic coordinates employed in this study have recently been published by Hsu et al. (PDB code 5UR9, 2.2 Å resolution)[19] and contained the FABP5 isoform cocrystallized with compound (S)-SBFI-26 originally identified by Berger et al.[12] Having a cocomplex is valuable from a structure-based screening standpoint because the target binding site is “preformed” to accommodate small molecule binding. Ligand (S)-SBFI-26 (hereafter referred to simply as SBFI-26) provides a key “reference” that can be used to identify additional compounds that make similar interactions with specific protein residues. Chain G protein from the oligomeric X-ray structure was used as the structure for docking because the ligand was well organized in the binding pocket and there were no missing residues or other factors that might negatively impact the virtual screen.

Docking Setup

Following previously reported procedures and protocols,[23−25] the FABP5-SBFI-26 cocrystal complex 5UR9[19] was prepared for docking. Briefly, coordinates were extracted from chain G of the structure, and the AMBER16[26] program antechamber and tleap modules were used to add hydrogen atoms to the protein, assign ff99SB[27] parameters to the protein, and assign GAFF[28] augmented by AM1-BCC[29,30] parameters for the ligand SBFI-26. Default AMBER protonation states were employed (Asp and Glu deprotonated, Lys and Arg protonated, His protonated at epsilon nitrogen). Hydrogen atom orientations were optimized using the sander module in AMBER16 for a maximum of 100 cycles of minimization, with heavy restraints (1000.0 kcal mol–1 Å–2) on all non-hydrogen atoms. Protein and ligand coordinates were extracted and saved in MOL2 format required by DOCK6. Finally, the reference ligand SBFI-26 was minimized in the context of the DOCK setup with positional restraints of 5 kcal mol–1 Å–2, resulting in a small 0.81 Å root-mean-square deviation (RMSD) from the crystal geometry (heavy atoms).

The FABP5 structure was then used to create a molecular surface using the program dms(31) which was used as input for the program sphgen(32) to generate docking spheres which will guide the initial orientation of ligand “anchors”, the largest rigid fragment from which the rest of a molecule can be rebuilt. Spheres within 8 Å of the minimized reference ligand (SBFI-26) were selected, up to a maximum of 75, which were used to define a bounding box (binding site) for docking by including an 8 Å margin in each dimension from the entire sphere group. The DOCK grid(33) module was then used to prestore van der Waals (VDW) and electrostatic (ES) terms for the protein at grid point locations which were 0.3 Å apart which speeds up the calculations. Following lab protocols, a 6–9 Lennard–Jones potential was used for the VDW term to soften the energy landscape and a distance dependent dielectric (ddd) = 4r was used for the ES term as a crude approximation of solvation effects. Figure visualizes the key components of the docking setup for FABP5.

Figure 1

Docking setup for virtual screening targeting FABP5 (PDB 5UR9),[19] protein surface in tan, docking spheres represented by blue spheres, SBFI-26 reference ligand in cyan, and docking bounding box in gray.

Virtual Screening Protocol

A library of ∼2 M druglike small molecules from the ZINC database[22] was docked (FLX protocol)[25] to the FABP5 active site using the program DOCK6.[23] The FLX protocol samples ligand torsion and rigid body degrees of freedom using the anchor and grow algorithm with the receptor (FABP5) being held rigid. We have previously employed this DOCK procedure to successfully identify inhibitors targeting HIVgp41,[34−37] FABP,[12,19] HER2,[38] and BoNT.[39,40] For each screened molecule, only the best grid-based energy conformation (pose) was retained. Each pose was subsequently minimized in DOCK Cartesian energy (DCE) space, using a 6–12 Lennard–Jones potential and a distance dependent dielectric constant (ddd = 4r), to remove grid-energy approximations and facilitate prioritization by enhanced scoring functions.

Enhanced Scoring Functions

The energy-minimized poses were evaluated and prioritized by a series of enhanced scoring functions available in DOCK6.8, including footprint similarity score (FPS)[41] score, pharmacophore matching similarity (FMS)[42] score, Hungarian matching similarity (HMS)[43] score, and volume overlap similarity (VOS) score. Footprints are a breakdown of the energetic interactions between a ligand and protein by primary sequence, and an FPS score between any docked candidate and a reference (in this case ligand SBFI-26 from the X-ray structure) is quantified by computing the Euclidean distance between the two interaction vectors (footprint interaction patterns). The FMS score quantifies similarity between docked candidates and a reference based on the number of matched pharmacophore features and quality of the matches. The HMS score uses the Hungarian algorithm[44,45] to calculate an RMSD-like metric which quantifies the overall minimum distance between two molecules based on comparisons of atom pairs of the same atom type. Lastly, the VOS score, based on the algorithm reported by Sastry et al.,[46] quantifies geometric volume overlap between candidates and a reference molecule using all-atom, hydrophobic and hydrophilic atoms, and positively and negatively charged atom definitions.

Compound Selection (Initial Screen)

To promote molecular diversity among the compounds selected for purchase, the top 100 000 molecules ranked initially by DCE score were clustered using a best first clustering algorithm with a Tanimoto similarity score cutoff of 0.95 based on the MACCS[47] fingerprint method as implemented in the Molecular Operating Environment (MOE) software suite.[48] To help prioritize compounds, the resultant clusterheads were then reranked by eight different criteria: (1) DCEVDW+ES, DOCK Cartesian energy score consisting of van der Waals plus electrostatic terms, (2) FPSVDW+ES, footprint similarity score consisting of van der Waals plus electrostatic terms, (3) FPSVDW, footprint similarity score consisting of only the van der Waals term, (4) FPSES, footprint similarity score consisting of only the electrostatic term, (5) Total Score, a linear combination of DCEVDW+ES and FPSVDW+ES, (6) FMS, pharmacophore matching similarity score, (7) HMS, Hungarian matching similarity score, and (8) VOS, volume overlap similarity score. Additional molecular descriptors including number of rotatable bonds, molecular weight, numbers of hydrogen bond donors and acceptors, number of chiral centers, SlogP, formal charge, and logS were calculated (MOE and/or DOCK programs) and used in some cases to eliminate compounds with undesirable properties. Figure outlines the overall virtual screening protocol. Upon visual examination of docked poses in the binding site from different rank-ordered lists, 78 compounds were ultimately purchased (ChemDiv vendor) for experimental testing.

Figure 2

DOCK6 Virtual Screening Protocol

Fluorescence Displacement Binding Assays

Purification of recombinant human FABP3 and FABP5 were performed as previously described.[11] Binding assays were carried out in 96-well Costar plates (Corning Life Science, Kennebunk, ME). Recombinant human FABP3 or FABP5 (3 μM) was incubated with the fluorescent probe (500 nM) in binding assay buffer (30 mM Tris-HCl, 100 mM NaCl, pH 7.6). Competitor test compounds (0.1–50 μM) were then introduced to the well and mixed, and the system was allowed to equilibrate for 20 min at 25 °C in the dark. All experimental conditions were tested in triplicate. Each independent assay also included a strong competitive binder (arachidonic acid, 10 μM) as a positive control for probe displacement and background wells that did not contain any protein. Loss of fluorescence intensity was monitored with an F5 Filtermax Multi-Mode Microplate Reader (Molecular Devices, Sunnyvale, CA) using excitation (ex.) and emission (em.) wavelengths appropriate for each probe (NBD-stearate ex./em. = 465/535 nm, DAUDA ex./em. = 345/535 nm). Single point experiments utilized the NBD probe while dose response experiments utilized the DAUDA probe. Following background subtraction, raw fluorescence intensity values were normalized and fit to a one-site binding analysis using the Graphpad Prism software (Prism version 7.0 for Mac OS; Graphpad Software Incorporated, La Jolla, CA).

Cytotoxicity

Human umbilical vein endothelial cells (HUVEC) were grown in Endothelial Cell Growth Medium (Sigma-Aldrich) supplemented with 10% FBS, 100 units/mL penicillin/streptomycin, and 1 mM sodium pyruvate. Cell viability was assessed by MTT colorimetric assay. HUVEC cells were seeded into 96-well culture plates and grown in complete medium until the following day. The media was aspirated from the wells, gently washed with PBS, and 200 μL serum-free medium supplemented with 0.05% bovine serum albumin (BSA) and containing indicated concentrations of compound or vehicle (0.5% DMSO). Each plate contained three wells with 0.5% sodium dodecyl sulfate (SDS) that was used as a positive control for cell death and three wells with no cells seeded to be used as a background reading. The plates were then incubated for 24 h at 37 °C, at which time the drug-containing media was removed and 200 μL serum-free media containing 0.5 mg/mL MTT was added to all wells and incubated for 3.5 h at 37 °C. The media was gently aspirated and DMSO was added to solubilize the resulting formazan crystals. Plate absorbance was read at 585 nm with an F5 Filtermax Multi-Mode Microplate Reader. Background readings were subtracted and the IC50 was determined by plotting the resulting curve with Graphpad Prism software.

Analog Selection (Secondary Screen)

To further interrogate the scaffolds of experimentally verified active compounds, hits from the first screen were used to select analogues via two different methods: (1) The first method involved rescoring the originally docked library with FPS, FMS, HMS, and VOS similarity-based functions using the most promising hit compound as a reference. Compounds related to the hit were chosen from the intersection of the different rank-ordered lists (100 top-scoring molecules each). (2) The second method employed the similarity search feature in ZINC15 to identify commercially available analogues having a Tanimoto score of 0.75 or higher. Available compounds were downloaded, and flexible docking and continuous space minimization were performed in a similar manner as the original virtual screen. To help generate a “consensus pose” the DOCK multigrid scoring function[49] was used which incorporates a footprint-based energy term which favors geometries making similar interactions as a reference (discussed in Results and Discussion).

Molecular Dynamics Simulations and Free Energy Calculations

Molecular dynamics simulations and free energies of binding (MM-GBSA method[50,51]) were also performed to further interrogate the geometric and energetic stability of the experimentally verified hits. Starting from each DOCK-predicted pose, MD-ready complexes were constructed using the AMBER16 package.[26] Briefly, proteins were parametrized with the ff14SB force field,[52] and ligands were parametrized using GAFF[28] (augmented with AM1-BCC[29,30] partial charges) which were assigned using the antechamber module. Each complex was solvated with TIP3P[53] water in an octahedron with a 13 Å margin in each direction. For all neutral ligands, one sodium counterion was added to the systems to keep the formal charge of the system neutral. No additional salts were added.

A previously employed nine-step minimization and equilibration protocol[37] was used to relax each system prior to production MD. The first steps employed energy minimization of hydrogens atoms and solvent molecules with solute heavy atoms heavily restrained (20 kcal mol–1 Å–2) to their initial starting positions, followed by a minimization of the entire complex. Each system was then heated sequentially from 50 to 300 K, over 250 ps, with the same 20 kcal mol–1 Å–2 solute restraints. Five additional MD equilibrations were performed (4 × 200 ps, 1 × 500 ps) with positional restraints on protein backbone and ligand atoms gradually reduced from 5.0 kcal mol–1 Å–2 to 0.1 kcal mol–1 Å–2 with the last equilibration step having no ligand restraints. Lastly, 20 ns of production data was collected for each system using the CUDA-accelerated version of pmemd in AMBER16 under NPT conditions with a 0.1 kcal mol–1 Å–2 positional restraint on the protein backbone only. Snapshots of each trajectory were saved every 5 ps for data analysis and fit to the initial MD frame using protein backbone heavy atoms. We have employed similar MD protocols (i.e., overall setup, equilibration and production procedures, simulation lengths) to characterize protein–ligand binding in a variety of systems.[35−38,40]

Analyses of the trajectories including root-mean-square deviation (RMSD) calculations, distance measurements, and clustering which were performed with the cpptraj(54) module. For RMSD calculations for ligands no additional fitting was performed. Clustering of evenly spaced snapshots from trajectories were performed using the hierarchical agglomerative algorithm as implemented in cpptraj. End-state free energies of binding (ΔGbind) were estimated using the “single-trajectory” MM-GBSA[50,51] method, facilitated by the AMBER16 program MMPBSA.py,[55] using the periodically saved snapshots. Error analysis was performed using autocorrelation functions (ACF) and block-averaged standard errors of the mean (BASEM)[56,57] as previously reported by our group[58,59] (see the Supporting Information). The MD snapshots were also used to compute time-averaged molecular footprints as outlined in previous work.[37,40,59−61]

Results and Discussion

Virtual screening outcomes: Overlap between highly ranked compounds

As described in Methods, compound prioritization involved consideration of eight different rank-ordered lists of compounds generated using different DOCK scoring functions that include similarity-based metrics (Figure . Numerous docked candidates might have been expected to be highly ranked by different functions and therefore would be included in multiple lists, given that only a single reference ligand was used (SBFI-26). However, since different scoring functions quantify 3D similarity from different perspectives (i.e., footprint similarity, pharmacophore, Hungarian, volume overlap), a significant level of diversity is in fact observed as shown in Table .

Table 1

Number of Top-Ranked Compounds in Common for Different Ligand Ensembles (N = 200 Each) Obtained by Different Scoring Functionsa

	DCEVDW+ES	FPSVDW+ES	FPSVDW	FPSES	TotalScore	HMS	FMS	VOS
DCEVDW+ES	200	0	0	0	61	0	0	0
FPSVDW+ES	0	200	37	55	5	3	11	8
FPSVDW	0	37	200	4	0	6	5	8
FPSES	0	55	4	200	6	1	11	3
TotalScore	61	5	0	6	200	0	1	0
HMS	0	3	6	1	0	200	12	19
FMS	0	11	5	11	1	12	200	13
VOS	0	8	8	3	0	19	13	200

DCEVDW+ES = DOCK Cartesian energy score consisting of van der Waals plus electrostatic terms, FPSVDW+ES = footprint similarity score consisting of van der Waals plus electrostatic terms, FPSVDW = footprint similarity score consisting of only the van der Waals term, FPSES = footprint similarity score consisting of only the electrostatic term, TotalScore = a linear combination of DCEVDW+ES and FPSVDW+ES, FMS = pharmacophore matching similarity score, HMS = Hungarian matching similarity score, VOS = volume overlap score. Matrix values quantify the number of top-ranked compounds in common for pairs of ligand ensembles (N = 200 each) obtained by two different DOCK scoring functions.

For example, among the top-200 ranked compounds, the DCEVDW+ES list (standard DOCK energy function) shows overlap only with TotalScore (N = 61) which is reasonable because DCEVDW+ES is included as component of TotalScore and it is an energy-based function only (not similarity-based). Although there is overlap between FPSVDW+ES with FPSVDW (N = 37) or FPSES (N = and 55), again because FPSVDW and FPSES both contribute to FPSVDW+ES score, there is little overlap between the components themselves (N = 4) because they quantify similarity only in terms of their van der Waals or electrostatic interactions patterns, respectively. The generally small overlap between the four main different similarity-based list (FPSVDW+ES, HMS, FMS, and VOS) indicates significant diversity among the different group compounds which highlights the potential benefit of using multiple scoring metrics to help prioritize compounds for purchase and experimental testing. On the other hand, the fact that there is some overlap between the different groups suggests there are a small number of compounds that were predicted to bind to FABP5 in an extremely similar fashion as the SBFI-26 reference with respect to multiple binding descriptors.

To visually emphasize how use of different functions leads to different outcomes Figure shows docked ensembles (200 clusterheads each) rank-ordered by (a) DCEVDW+ES, (b) FPSVDW+ES, (c) TotalScore, (d) HMS, (e) FMS, and (f) VOS. As observed in prior studies,[34,35,40] use of the DCEVDW+ES function (Figure a) generally leads to compounds that are larger in size (MW bias) compared to other ensembles. In contrast, the ligands prioritized using similarity-based methods are more compact (especially Figure b,d,f), and are more spatially clustered around the reference ligand SBFI-26 (not shown in figure for clarity). The TotalScore function was designed to provide more of a balance between a purely energetic (DCEVDW+ES) and a similarity-based method (FPSVDW+ES). The resultant TotalScore ensemble shown in Figure c does, in this case, appear to contain molecules that are more spatially balanced between those selected by either component alone (Figure a,b).

Figure 3

Overlay of top-ranked clusterheads (200 compounds each) prioritized by different scoring functions with FABP5: (a) DCEVDW+ES, (b) FPSVDW+ES, (c) TotalScore, (d) HMS, (e) FMS, and (f) VOS. See text for scoring function descriptions.

Virtual Screening Outcomes: Compounds Selected for Experimental Testing

To arrive at the prioritized group for purchase and experimental testing, the top-ranked clusterheads were visually inspected using 3D stereographics in the context of their predicted binding pose along with an examination of key molecular properties (descriptors) including: (1) position relative to the reference ligand, (2) similar footprint interaction pattern compared to the reference, (3) low number of chiral centers, (4) low computed ClogP values, and (5) Lipinski rules violations. In total, 134 compounds were selected for experimental testing divided roughly evenly among the different functions employed: DCEVDW+ES = 20, FPSVDW+ES = 28 (FPSVDW and FPSES were not explicitly used for compound prioritization due to their high overlap with the FPSVDW+ES ranked list, see Table ), TotalScore = 24, HMS = 26, FMS = 21, and VOS = 22. Of these, only 78 out of 134 were in stock and were purchased for experimental testing. One particularly notable feature of the FABP5 active site is that ligand binding (i.e., native substrate AEA or inhibitor SBFI-26[19]) involves a key electrostatic interaction with Arg129. Visualization of the 78 purchased compounds (Figure ) highlights this interaction in which different functional groups H-bond with Arg129 through carbonyl oxygens, nitrogens, sulfones, and carboxylic acids.

Figure 4

Close-up view of FABP5 residue Arg129 with the 78 docked and purchased compounds (orange) highlighting a key conserved H-bond interaction (magenta lines).

Experimental Testing of Compounds Based on the DOCK Virtual Screen

The 78 purchased compounds were tested for FABP5 and FABP3 affinity at two concentrations (20 μM and 5 μM). The activity trends were consistent for the compounds at both concentrations thus for simplicity only the results at 5 μM are discussed below. As shown in Figure a, four compounds displayed 60% or greater probe displacement at 5 μM with FABP5, which is a comparable range of relative binding affinity to the positive control SBFI-26 (Figure , green). However, the two compounds with the strongest affinity for FABP5 (V015-5448, F255-0055) also showed strong affinity for FABP3 and therefore were not pursued. The fourth strongest hit (colored orange in Figure , compound Y070-2541 (ZINC09463091), was ultimately selected for additional study based on its reasonable affinity for FABP5 and low affinity to FABP3 which suggests a favorable selectivity profile. The docked pose for Y070-2541 also showed high overlap in terms of favorable FPS (5.67) and HMS (−0.21) similarity scores relative to the inhibitor SBFI-26 reference.

Figure 5

Fluorescence displacement binding assay results for the first set of 78 compounds (5 μM) against (a) FABP5 and (b) FABP3. The results were arranged in the order of their mean % Fluorescence value against FABP5. Control compounds are arachidonic acid (red) and SBFI-26 (green).

Further investigation of Y070-2541 (hereafter called STK-0) at different concentrations confirmed it to be a strong binder to FABP5. As shown in Figure a, the Ki value of STK-0 with FABP5 was determined to be 5.53 ± 0.89 μM which is ∼6-fold less active than SBFI-26 (0.86 ± 0.18 μM) under the current assay conditions. Considering that SBFI-26 is one of the most potent FABP5 inhibitors previously reported and the fact that STK-0 has a novel chemotype, it is a reasonable starting point for further refinement to improve affinity. Figure b visually shows the significant overlap between predicted binding pose for STK-0 (orange) and the X-ray pose for SBFI-26 (green). Notably, the docked scaffold of STK-0 shares many of the important aromatic and polar features with SBFI-26. It is worth noting that STK-0 was ranked only 38 081 by the standard DCEVDW+ES score, thus it would never have been selected for experimental testing if only intermolecular energy (i.e., DCEVDW+ES score) was used as the criteria for prioritization. This dramatic observation highlights the utility of similarity-based scoring functions, in this case HMS[43] score, during compound selection.

Figure 6

Comparison of STK-0 and SBFI-26 showing (a) 2D structures, (b) activities with FABP5, and (c) binding site geometries (STK-0 orange, SBFI-26 green).

Identification and Selection of STK-0 Analogues

The most promising compound (STK-0) from the virtual screen was subsequently used as a reference to identify a second group of analogs for experimental testing using two distinct procedures (see Methods): (1) 3D search queries of the initially docked library using DOCK similarity metrics and (2) 2D search queries of ZINC using molecular Tanimoto similarity score.

For the 3D searches (Figure a), the existing docked library of ∼2 M compounds was reprioritized with four DOCK similarity-based functions and provides an example of data mining. Analogous to the methods used in the original virtual screen, the data mining employed FPS, FMS, HMS, and VOS scores; however, the reference was the docked pose of STK-0 (not SBFI-26). Ultimately, seven compounds were identified from the intersection of the four rank ordered lists (100 molecules each). As expected, visualization of this group (Figure a) showed good correspondence in terms of functional group overlap with the STK-0 reference.

Figure 7

Comparison of DOCK predicted binding geometries for parents (cyan) vs analogues (magenta) derived from (a) 3D similarity searches of the originally docked library (N = 7) and (b) 2D similarity searches in ZINC based on parent compound STK-0 (N = 48).

For the 2D searches (Figure b), STK-0 was employed to query the ZINC database, using a Tanimoto similarity cutoff of 0.75, which resulted in the identification of an additional related group of compounds (N = 88). It is worth noting that STK-0 and the identified analogues all have two chiral centers, thus among the 88 hits, there were actually only 22 chemically unique molecules (four enantiomers each). For any analogue to be considered further, we imposed the requirement that DOCK-generated poses for at least one of enantiomers be similar to that of the reference. This was deemed important, given that a congeneric series would be expected to have similar binding geometries to establish and further predict structure–activity relationships.

However, our initial examination of the 22 × 4 poses generated with DOCK using the standard DCE energy function yielded a wider diversity of geometries than expected, despite the molecules having a common core. We hypothesize this outcome was a result of the following: (1) the binding site of FABP5 is relatively large which, in some cases, allow flexible ligands to adopt a variety of conformations; (2) each of the 22 unique molecules had 4 enantiomers, some of which could not adopt a similar pose as the parent due to chirality; (3) the scaffold of STK-0 and its analogues is branched, and the branches share somewhat similar chemical functionalities, mostly aromatic rings. Therefore, different arrangements (i.e., alternative quasi-symmetric conformations) can result in similar scores when only the DCE function was used.

As an alternative protocol, we redocked the 22 × 4 analogues using the DOCK multigrid similarity (MGS) function which includes a footprint-based term to reward poses making similar interactions as a reference (in this case SBFI-26 from the X-ray structure). With this protocol, the consistency of predicted binding geometries increased significantly such that 48 out of the 88 top-scored poses resembled that of STK-0 by visual inspection (Figure b). Overall, the calculations showed that at least one enantiomer from each of the 22 analogues could bind to FABP5 in a manner similarly as the parent STK-0. In total, 26 out of the unique compounds identified from the combination of the secondary 3D and 2D similarity searches were available for purchase and ordered for experimental testing. This group of compounds was coded STK-1 thru STK-26.

Experimental Testing of Analogues Derived from the Similarity-Based Searches

The STK group was subsequently tested using the same assay with the NBD-stearate probe at both 5 μM and 20 μM concentrations. Compounds were first reconstituted in DMSO to a concentration of 5 mM. Out of the 26 compounds, 4 were not soluble at the stock concentration and were not considered further. Figure a shows results for 22 compounds at 5 μM in comparison to the parent compound STK-0 and the previously identified inhibitor SBFI-26. Encouragingly, two compounds (STK-15, STK-22) displayed better affinity for FABP5 than STK-0 and one compound (STK-15) appeared to be slightly more potent than SBFI-26. In terms of selectivity for FABP5 over FABP3 (Figure b), although no analogues were as selective as the STK-0 parent, STK-15 and STK-22 were more selective than SBFI-26.

Figure 8

Fluorescence displacement binding assay results for STK compounds (5 μM) against (a) FABP5 and (b) FABP3. The results were arranged in the order of their mean % Fluorescence value against FABP5. Control compounds are arachidonic acid (red) and SBFI-26 (green).

Structurally, all analogues share the same central dihydropyrrole group that is predicted to maintain the key interaction with Arg129 observed in the SBFI-26 X-ray structure. The analogues contain modifications to the STK-0 ethyl-phenyl, chloro-phenyl and benzofuran groups that include substitutions with functional groups of similar size, replacement of chlorine with other halogens or larger functional groups, and/or the addition of other functional groups to either phenyl ring. Remarkably, within the STK series rank-ordered by their %fluorescence with FABP5 (Figure a), the 14 top ranked compounds (STK-15 through STK-20) all retain the N-substituted ethyl-phenyl group relative to the STK-0 parent, with the sole exception of STK-16 (methyl-phenyl) which interestingly has the largest error among this group. Conversely, the eight compounds showing much weaker activity for FABP5 (STK-24 through STK-2) do not have an N-substituted ethyl-phenyl with the exception of STK-12 and STK-13. Taken together, these SAR trends strongly indicate that the N-substituted phenyl group should be conserved and the chloro-phenyl group is a promising position for further exploration. The fact that the most potent analogues (STK-15 and STK-22) replace the chlorine with oxy-phenyl or methoxy-phenyl functionality, two bulkier groups unique among the series of compounds tested, suggest that exploration of alternative bulkier aromatic rings would be worthwhile.

For the three STK compounds with the most FABP5 activity at 5 μM (STK-15, STK-21, and STK-22), rigorous dose response affinity and cytotoxicity experiments were subsequently performed as shown in Figure . Encouragingly, all three compounds showed clear dose–response behavior, had excellent cytotoxicity profiles (>75% cell viability up to 100 μM), and the measured Ki values for STK-15 (1.40 ± 0.35 μM), STK-21 (3.13 ± 0.21 μM), and STK-22 (2.90 ± 0.93 μM) were lower than the initial parent STK-0 (5.53 ± 0.89 μM).

Figure 9

2-D structures and dose response curves for binding (blue) and HUVEC cell viability (red) for (a) STK-15, (b) STK-21, and (c) STK-22. Values represent the average of three independent experiments.

The four hits were also examined with respect to possible nonspecific effects as a result of colloidal aggregation, pan-assay interference compound (PAINS) liabilities, or promiscuity. According to the Aggregator Advisor (advisor.bkslab.org) Web site, none of the hits had been previously reported as an aggregator. Conversely, the ZINC database (zinc15.docking.org) searches yielded 2-D similarity scores ranging from 0.55 to 0.60 to a previously reported aggregator ZINC13127469. However, the fact that the four hits, in particular STK-0 and STK-15, preferentially bind to FABP5 over FABP3 (Figure ) suggests their measured activities are not due to colloidal aggregation. The well-behaved dose–response curves in comparison to the known inhibitor SBFI-26 (Figure b), tested under the exact same conditions, also suggests that activity is not a result of aggregation. According to cbligand (cbligand.org/PAINS) and SwissADME (swissadme.ch), none of the four hits contain PAINS liabilities. Finally, in terms of promiscuity, PubChem (pubchem.ncbi.nlm.nih.gov) searches did not show that the compounds had previously been tested.

Ensemble-Based Characterization of STK Enantiomers

It should be noted that the hits share a common scaffold (Figure a, Figure a–c); however, as the experimentally tested samples were racemic mixtures (four enantiomers each) it is not known if different enantiomers might have similar activities. To explore the effects of stereochemistry on binding in greater detail, we performed molecular dynamics (MD) simulations and free energy calculations for each enantiomer of each hit (4 × 4 setups). For each setup, four MD replicas were executed following the protocols outlined in Methods. Examination of the initial DOCK poses for the most potent compound STK-15 showed that top scored poses for 3 out of the 4 enantiomers share an overall consistent binding pose with considerable overlap in terms of functional group placement (Figure a). Furthermore, and in agreement with the trend obtained previously for two different enantiomers of SBFI-26,[19] the footprint patterns between the S,R-STK15 and R,R-STK15 enantiomers studied here also show striking accord (Figure b). Also, for the second most active compound STK-22, 2 out of the 4 enantiomers had top scored poses that shared this same geometry. Although the predicted best poses for STK-0 and STK-21 showed more variability, with only one enantiomer in each case adopting a similar pose, a consensus-like geometry was always within the top 5 DOCK results with scores well within the margin of error compared to the lowest energy. Thus, all MD simulations were initiated using poses that resembled this overall consensus pose.

Figure 10

(a) Overlay of 3 of the 4 STK-15 enantiomers: S,R-STK-15 (blue), S,S-STK-15 (tan), and R,R-STK-15 (magenta). (b) Footprint (per-residue energy) comparison between S,R-STK-15 (blue) and R,R-STK-15 (magenta). Energies in kcal/mol.

Table summarizes the MD results which include (1) ligand heavy atom RMSD values to their respective initial geometries (RMSD to pose), (2) free energies of binding (MM-GBSA method), and (3) distances between the ligand dihydropyrrole scaffold and Arg109 (defined in Table legend, potential H-bonding illustrated in Figure . In general, the simulations showed good ligand stability with 12 of the 16 enantiomers yielding RMSD values under 3.0 Å to their initial DOCK consensus pose (values ranged from 1.67 to 3.84 Å). The computed free energies of binding were also favorable (i.e., negative) although there was a relatively wide range of values (−26.31 to −45.00 kcal/mol). The initial experimentally verified hit (S,R-STK-0, Table , bolded cells) yielded a 2.76 Å RMSD and a free energy of binding of −32.12 kcal/mol. For comparison, the known inhibitor SBFI-26, starting from its crystallographic binding geometry, had an RMSD of 1.87 Å ± 0.46 and a free energy of binding of −33.83 ± 0.54 kcal/mol, which suggests the specific MD protocols used here are robust. Overall, the simulation results reaffirm the predicted consensus binding geometry in Figure a.

Table 2

Molecular Dynamics Simulation Results for Four Enantiomers of STK-0, STK-15, STK-21, and STK-22a

	R,S-STK-0	R,R-STK-0	S,S-STK-0	S,R-STK-0b	averages
RMSD to pose (Å)c	1.67 ± 0.42	2.77 ± 1.28	3.20 ± 0.31	2.76 ± 0.40	2.60
ΔGbind (kcal/mol)d	–27.56 ± 0.25	–30.78 ± 0.54	–30.91 ± 0.16	–32.12 ± 0.40	–30.34
Lig–Arg distance (Å)e	5.16 ± 0.85	5.82 ± 1.87	7.15 ± 0.85	5.44 ± 0.67	5.89

The results for each enantiomer were averaged over four individual MD replicas of 20 ns each.

S,R-STK-0 (bolded cells) was the enantiomer first docked to FABP5 and the original hit.

Reference based on consensus pose. Fluctuations in standard deviation.

Single trajectory MM-GBSA method without entropy. Fluctuations in block-averaged standard error of the mean at a block size of 200 frames (see the Supporting Information).

Lig–Arg distance defined as the distance between the ligand dihydropyrrole scaffold oxygen adjacent to N and the closest nitrogen on Arg109. Fluctuations in standard deviation.

Figure 11

Comparison of poses for (a) R,R-STK-15 and (b) S,R-STK-15 from representative snapshots from the largest MD clusters observed over 80 ns of simulation.

Interestingly, the two most potent compounds STK15 and STK22 have, on average over all enantiomers (Table , far right column), more favorable predicted free energies of binding (−36.39 and −38.34 kcal/mol) compared to STK0 and STK21 (−30.34 and −29.49 kcal/mol). Also, for STK15 and STK22 in particular, the R,S (−40.57 and −42.35 kcal/mol) and R,R (−36.95 and −45.00 kcal/mol) enantiomers are predicted to bind much more tightly than their respective S,S (−31.77 and −32.62 kcal/mol) or S,R (−36.25 and −33.39 kcal/mol) forms. An inspection of the MD trajectories revealed an interesting pattern for the R forms (i.e., R,S, and R,R forms) of STK-15 and STK-22 in which the ligand would shift slightly in the binding pocket which enabled additional H-bonding with Arg109 (Figure . Compellingly, calculated ligand–Arg109 distances (Table ) for these four complexes were significantly smaller (2.91–3.84 Å) than the other 12 systems (5.16–7.15 Å). Thus, the more favorably calculated free energies of binding observed for these four enantiomers appear to be a direct consequence of increased H-bonding with Arg109. Together, the analysis suggests that enantiomerically pure compounds could show increased activity.

Origins of Enhanced H-Bonding with Arg109 for R-STK-15 and R-STK-22

To understand why only the R forms (R,S, and R,R) of STK15 and STK22 showed enhanced H-bonding with Arg109, we compared the initial starting coordinates with those from MD trajectories. Figure a illustrates that the initial consensus poses for R,R-STK-15 (meta O-phenyl, orange) and R,R-STK-0 (meta chlorine, gray) were well overlaid, particularly along the central pyrrole ring and that the ligand pyrrole oxygen was not within H-bonding distance of Arg109. All enantiomers had similar starting poses. Thus, as expected, there was no initial structural biases that might have favored H-bonding with Arg109 for certain ligands or enantiomers and any conformational shifts happened spontaneously during the MD simulations.

Figure 12

(a) Starting coordinates for R,R-STK-15 (meta O-phenyl, orange) and R,R-STK-0 (meta chlorine, gray) in the FABP5 X-ray crystal structure. (b) Comparison of starting coordinates (orange = R,R-STK-15, purple = Phe65) vs representative MD snapshot (green = R,R-STK-15, cyan = Phe65).

One standout feature observed during the MD trajectories was a relatively minor rotation of the FABP5 Phe65 side chain (cyan vs purple) as illustrated in Figure b for R,R-STK-15. Mechanistically, the Phe65 rotation appears to accommodate an extended conformation of the N-substituted ethyl-phenyl group (green) on R,R-STK-15 into a groove formed by Ile54, Thr56, Phe65, Gln96, and Arg109, which enables the ligand to shift to the left (green vs orange) and interact with Arg109. Concurrently, this shift appears to place the larger phenoxy-phenyl (STK-15) or phenoxy-ethylphenyl (STK-22) groups at an optimal distance to interact with a pocket formed by Met25, Val28, Leu32, and Lys61. Due to their smaller sizes, ligands STK-0 (chlorophenyl) and STK-21 (bromophenyl) were not able to reach both this pocket and Arg109 simultaneously. In addition, as the initial X-ray conformer of Phe65 (purple) would clash with the N-substituted ethyl-phenyl group conformation seen if H-bonding with Arg109, this explains the absence of DOCK poses that interact with Arg109. Enhanced H-bonding and binding energy seen in MD simulations of R,S-STK-15, R,R-STK-22, R,S-STK-22 are expected to involve the same proposed mechanism.

A natural question is why do the S forms (i.e., S,S or S,R) not shift to interact with Arg109 (distances all greater than 5 Å)? Figure a compares representative MD snapshots between S,R-STK-15 (orange) and R,R-STK-15 (green) in which the pyrrole scaffold of R (green) shifts slightly (∼1.2 Å) to the left of S (orange) to facilitate H-bonding with Arg109. In contrast, as a result of the flipped stereochemistry of the phenoxy-phenyl group in the S form, the N-substitute phenyl group is now positioned upward relative to the pyrrole ring which would clash in the protein conformation adapted to R (Figure b, black arrow, blue surface clash). On the other hand, the reverse experiment of placing R,R-STK-15 in the protein conformation coupled to the S enantiomer shows no such clash (Figure c) and ligand maintains most of the same interactions with the site. Taken together, the differences inherent to the 3D spatial arrangement of functionality in R vs S forms, coupled with the extra stabilization available as a result of larger functionality in STK-15 and STK-22 (phenoxy-phenyl or phenoxy-ethylphenyl) compared to STK-0 and STK-21 (Cl or Br), enables enhanced H-bonding with Arg109 when Phe65 undergoes a rotameric shift.

Figure 13

(a) Binding site comparison between S,R-STK-15 (orange ligand, tan protein) and R,R-STK-15 (green ligand, blue protein). (b) Overlay showing the S,R-STK-15 ligand (orange) in the R,R-STK-15 protein (blue) conformation through protein backbone alignment. (c) Overlay showing the R,R-STK-15 ligand (orange) in the S,R-STK-15 protein (blue) conformation.

Conclusion

The primary goal of this work was to identify small organic molecules with inhibitory activity and specificity for FABP5 compared to FABP3. Employing a structure-based screening strategy (DOCK6 program), ∼2 M commercially available compounds were docked to the FABP5 active site structure originally cocomplexed with SBFI-26 (a compound previously identified by our group through virtual screening and confirmed by X-ray crystallography (Figure ). The 100 000 top-scoring compounds (DCE function) were retained and used to create additional smaller independent rank-ordered lists for compound prioritization based on different scoring metrics that employed unique similarity-based functions (Figure ). There were no overlapping results between the top 200 compounds prioritized by energy-based (DCE) vs similarity-based methods (FPS, HMS, FMS, VOS), and there were relatively few overlaps using different similarity-based scoring metrics (Table ). This helps validate the hypothesis that use of multiple scoring metrics leads to a more diversified pool of candidates and that the different similarity-based functions encode different information. Ultimately, 78 compounds from the initial virtual screen were purchased for experimental testing (Figures and 5).

Using a fluorescence displacement binding assay with an NBD-stearate substrate, 4 out of the 78 tested compounds showed better than 60% fluorescence at 5 μM concentration against FABP5. One compounds in particular (Y070-2541, code name STK-0), identified through Hungarian Matching Similarity (HMS) scoring, had a good selectivity profile against FABP3 (Figure ) and a reasonable Ki value (5.53 ± 0.89 μM) compared to the known control SBFI-26 (0.86 ± 0.18 μM) (Figure ) and was thus retained for further analysis. Subsequent similarity searches using STK-0 as reference, with the additional requirement that the analogues should have the same general binding pose as the parent (Figure ), yielded 26 additional compounds that were purchased (STK-1 to STK-26 series). All of the analogues contained the same central dihydropyrrole scaffold as the STK-0 parent. Encouragingly, several of the analogues showed better or comparable levels of activity to STK-0 (Figure ), which helps validate our hypothesis that chemically similar compounds would yield comparable biological activities.

The results also allowed some preliminary SAR to be determined including the observation that the N-substituted phenyl group should likely be conserved and that the meta chloro-phenyl group is a promising position for further exploration. Of the four most active hits (STK-0, STK-15, STK-21, and STK-22), the two most potent analogues (STK-15 and STK-22) replaced the chlorine with bulkier oxy-phenyl or methoxy-phenyl functionality suggesting that further explorations at this position would be worthwhile. Notably, the most potent hit STK-15 (1.40 μM) was only 2-fold less active than the SBFI-26 control (0.86 μM).

Somewhat surprisingly, although the DOCK6 consensus binding geometry for the STK series was similar to the crystallographic pose of SBFI-26, follow up MD simulations revealed that only the R form (chiral center on central scaffold) for two of the four hits (STK-15 and STK-22) consistently underwent a subtle but important conformational change that enabled new H-bonds to be formed with Arg109 (Table , Figure ). This new interaction was in addition to the canonical Arg129 H-bond observed in X-ray structures with fatty acid substrates or SBFI-26. Subsequent computational analysis led to the conclusion that only the R forms of STK-15 and STK-22 interact with Arg109 given inherent geometric differences arising from R vs S stereochemistry, and their larger size compared to STK-0 and STK-21 which stabilizes binding following a rotameric change in the side chain of Phe65 (Figures –13). The analysis suggests that enantiomerically pure R,R-STK-15 and R,S-STK-15 may show enhanced activity against FABP5, and this aspect would be worthwhile to investigate in future work. In summary, this study has demonstrated the ability of similarity-based virtual screening methods to identify experimentally validated compounds against FABP5 and provides a strong starting point for further optimization efforts to improve their activity.