A Computational Approach for Identifying Experimental or Approved Drugs That Can Be Repurposed for the Treatment of Type-2 Diabetes

Abstract

Approved and experimental drugs can be utilized for new indications as illustrated in the case study presented herein. In this case study, allopurinol (trade name Zyloprim and Aloprim) which is currently utilized for gout, was retrieved from the Drug Bank and evaluated for a new indication. Utilizing a catechin derivative as a scaffold, a derivative was designed incorporating allopurinol. This novel molecule was predicted to act as an allosteric inhibitor of fructose 1,6-bisphosphatase (FBPase), a control point for entry into the biochemical pathway gluconeogenesis. The predicted inhibition was validated with a colorimetric assay. Potential toxicity was assessed using a HepG2 MTT assay. As an inhibitor of this enzyme, the novel molecule proved to be both potent and non-toxic in cell-based assays. Once optimized and tested in vivo, the novel molecule may be potentially used as a therapeutic agent for type-2 diabetes mellitus inhibiting FBPase. This action prevents the de novo synthesis of glucose and potentially contributes to lowering blood glucose levels for patient populations that are genetically prone to chronic high blood glucose leading to insulin resistance. The computational approach to the design of the novel potential lead compound is discussed in detail and validation data presented.

1. Introduction

Type-2 diabetes mellitus is recognized as a global epidemic by the World Health Organization with an estimated 500 million affected worldwide. The Center for Disease Control (CDC) estimates that 30 million in the US are affected. Most of this is caused by lifestyle and improper diet, a direct result of the mismanagement of agricultural and food distribution systems. Certain populations, such as Native Americans and Polynesians, appear to be predisposed to type-2 diabetes. This might be a result of genetic adaptations brought on by chronic food scarcity [1]. In both cases, type-2 diabetes is characterized by impaired insulin sensitivity or availability and increased endogenous glucose production (EGP). The primary source of EGP is gluconeogenesis in the liver, which is typically three-fold greater in type-2 diabetics. It is widely recognized that when gluconeogenesis is curtailed, this provides a valuable therapy for type-2 diabetes [2]. Analysis of the gluconeogenesis pathway has suggested that the best target for inhibition is the fructose 1,6-bisphosphatase (FBPase) enzyme. Yet after many years of work, there is no FBPase inhibitor that has reached the market. The FBPase molecule is a homotetramer, each monomer with an active site and two allosteric binding sites. The active site contains highly conserved amino acids and so is not a suitable drug target. The AMP binding site is an allosteric site, which has been known and targeted for many years for drug development [1, 2, 3, 4]. A novel allosteric site was discovered at the tetramer interface by Pfizer [5]. Although a potent inhibitory site, the volume of the pocket is made up of contributions from the four monomers and is highly mobile, so there is difficulty in predicting the actual binding modes of selected molecules. In addition, a 2:1 molar ratio of drug to protein is often required for effective inhibition, which makes drug studies challenging. This study focuses on the FBPase protein target AMP binding site relying heavily on similarity searches with the natural inhibitor AMP.

Herein is a case study focused on the early stages of the development of an FBPase inhibitor (lead compound) that targets the adenosine monophosphate (AMP) allosteric binding site is presented. This has led to the computational identification and assay verification of several effective compounds in the Drug Bank database. The compounds are currently approved drugs. If any of these drugs can be fully verified to be effective for a new application, this could lead to a breakthrough treatment for type-2 diabetes mellitus.

2. Strategies and methods

As shown in Figure 1, the repurposing workflow is a linear process. The process begins with literature searches then similarity searches of known inhibitors with searches of the Drug Bank database [6]. Once compounds are identified, molecular docking is done using AutoDock Vina for both fixed and flexible (side chain) docking protocols [7]. Molecular dynamics follows for accurate physiological check on the binding score of the ligand/protein complex [8] followed by laboratory validation [9]. In the hit-to-lead optimization stage, there was a reiterative process of modifying the hit repurposed compound (treated as a scaffolding compound) with substituted functional groups. This process refines the molecule eliminating derivatives issues such as low solubility, cellular toxicity, and PAINS (pan assay interference screening) before proceeding to optimize the lead compound for in vivo animal studies and finally clinical trials.

Figure 1.

A virtual drug development platform (VDDP): a schematic representation of repurposing workflow.

Key research performed already published that supports this strategy includes:

1. The ZINC15 database of 3D structures of small molecules [10].

2. The Drug Bank database of approved and experimental drugs [6].

3. The 3D structures of human and pig liver fructose 1,6-bisphosphatase (FBPase), PDB code 1FTA [11], and 1KZ8 [5] for FBPase in complex with AMP.

4. Colorimetric assays able to confirm inhibition of FBPase [12, 13].

5. Site-directed mutagenesis studies on FBPase used to validate binding site residues and significant residues for allosteric mechanism [11, 13, 14].

6. The recognition of the AMP binding site on fructose 1,6-biphosphatase as an important allosteric drug target for type-2 diabetes [2, 3, 4, 11, 12, 13, 14].

7. A 2- stage process of docking and molecular dynamics to screen millions of compounds [15].

Research conducted by symmetric computing supports this case study includes the development of the virtual drug discovery platform (VDDP), which performs high throughput virtual screening by combining molecular docking with molecular dynamics. Screening of millions of small molecules from the ZINC15 database was done using VDDP, in order to identify high-scoring binders to the AMP binding site on fructose 1,6-biphosphatase. Identification of high-scoring binders that are Phase I approved deployed or experimental drugs that were also included in the Drug Bank database were cross-referenced.

3. Flexible docking and MD protocols

3.1 Preparation of small molecules and protein target for MD and geometric docking

Once repurposing molecules for the target were identified through literature searches to be the appropriate for alternative medicinal activities, the PubChem and Drug Bank database was searched these compounds and derivatives. Following identification, the total polar surface area (TPSA) was calculated with SwissADME [16]. Small molecules with druggable TPSA of approximately 90–120 A^2 [16, 17, 18] were fed into the SwissADME program, followed by the pkCSM program to screen for central ADMET properties to test for predicted toxicity using graph-based patterns [17]. The protein target binding site the FBPase AMP binding site was then selected based on geometric attributes and the docking score of the Drug Bank scaffolds. More specifically, the geometric attributes of the protein binding site coordinates were selected according to the ratio of molecular surface volume of the cavity over the total molecular surface volume of the protein. Ratios below 0.15 A^3 were not considered druggable target binding sites [19]. In addition, the docking score of the ligand to the protein was considered reasonable when less than ~50 micromolar. Once the protein target binding site was correctly identified, docking was performed with all of the selected drugs for repurposing (n = ~525) with both flexible and rigid residues. After initial docking (rigid) followed by flexible docking, a molecular dynamics protocol in NAMD was executed to check contacts and delta G of binding scores. Since convergence with flexible docking scores and MD scores occurred, flexible docking was used for drugs that were then altered with specific functional groups for further protein/drug docking studies.

The initial workflow for ligand and FBPase protein target preparation is shown in Figure 2. Two protein target PDBs that were chosen of FBPase were 1KZ8 [5] and 1FTA [2] to capture the AMP binding site and the dimer and tetramer interfaces in different starting conformations.

Figure 2.

Ligand/protein target computational platform workflow chart.

A guide of the initial identification of potentially druggable binding sites for a given target and previously approved drugs appropriate for drug repurposing studies.

Initially, the RCSB was searched for human and pig kidney structures of FBPase as these species have >85% identity and ~95% homology of amino acid sequences. Although the location of the AMP binding site of FBPase is known, the conformation of the amino acid side chains that make up the “gateway” and interior of the binding site change position when substrate or product is bound in the active site. The enzyme has two canonical states the T state known as inactive state and the R state the active state. The enzyme rotates around the dimer interface of the tetramer. Since the AMP binding site is on exterior side of the homotetramer on each monomer individual amino acid side chains shift position during this rotation.

The repurposed drugs were selected based on similarity searches to natural allosteric inhibitor in both PubChem and Drug Bank databases. Once the protein target center of position coordinates are selected and the repurposed drugs identified, static docking followed by flexible docking was done with the FBPase-repurposed drug complex. After molecular docking, the coordinates were run in MD protocols to check overall energy of the system and calculate the binding energy of the protein-drug complex. From the output from the VDDP, repurposed drugs from the Drug Bank were chosen for laboratory validation, derivatized as needed to enhance binding and the process of optimization towards lead compounds was done. In addition, the processed target protein can be “recycled” in its altered conformation out of MD to be tested with other repurposed drugs. In a reiterative process tailored to each small molecule, different amino acids can be made flexible prior to computationally expensive MD, simulating an induced fit model of protein-ligand binding using the platform. The reiterative process of flexible docking coupled with MD protocols decreases the computational time to generate new protein conformers significantly.

4. Molecular dynamics

4.1 MD simulations

Preparations for MD simulations were done using SwissDock parameter files for the catechin derivatives for MD [20]. The protein structure files were generated with docked ligands in VMD [21]. Molecular dynamics simulations were performed with the NAMD Program [8] using CHARMM37 parameters along with complementary CHARMM General Force Field optimized for ligand parameter files for each protein target/catechin derivative complex [22]. Periodic boundary conditions were imposed and explicit solvent/implicit solvent boundaries defined using a truncated rectangular box, ensuring a solvent shell of at least 10 A around the solute surrounding the protein-ligand complex. The solute was neutralized with potassium ions (K⁺/Cl⁻ ion pairs) to a concentration of 150 mM. The ions were initially placed at random, but at least 5 A from ligands and 3.5 A from one another. The resulting systems contained between 10,500 and 15,250 water molecules, corresponding to a total of 30–45,000 atoms.

Simulations employed periodic boundary conditions and electrostatic interactions were treated using the particle-mesh Ewald algorithm [23, 24] with a real space cutoff of pair list was built with a buffer region, and a list update was triggered whenever a particle moved by more than 0.5 A with respect to the previous update. Each system was initially subjected to energy minimization with harmonic restraints of 2 kcal mol⁻¹ A ̊⁻² on the solute atoms. The system was then heated to 310 K at constant volume during 100 ps. Constraints were then relaxed from 5 to 1 kcal mol⁻¹ A ̊ ⁻² during a series of 1000 steps of energy minimization (500 steps of steepest descent and 500 steps of conjugate gradient) followed by 50 ps of equilibration with restraints of 0.5 mol⁻¹ A ̊⁻² and 50 ps without solute restraints. The 50 ns production simulations were carried out at constant temperature (300 K) and pressure (1 bar) with a 2 fs time step. This was reiterated 10×. During these simulations, pressure and temperature were maintained using the Berendsen algorithm [25] with a coupling constant of 5 ps, and SHAKE constraints were applied to all bonds involving hydrogens [26]. Conformational snapshots were saved for further analysis every 10 ps. For comparison purposes, the isolated catechin complexes from each complex were also simulated alone using an identical protocol, creating a second set of ten 50 ns trajectories.

4.2 Potential energy of binding and free energy of binding calculations for MD output

The overall potential energy of protein-ligand system was evaluated by adding the final potential energy of the small molecule alone, and the final potential energy of the protein alone run through the same MD protocol. These values were compared with the protein/ligand complex. An expected, drop in the potential energy of the complex was confirmed. Absolute binding free energy (ABFE) calculation, with total annihilation of the ligand in the binding pocket followed by its reappearance at bulk state where the target protein was absent, was part of this process. As free energy is a state function, the alchemical FEP route to getting binding free energy of these derivates was as follows: (1) ‘locking the ligand’, restraining conformational, translational, and rotational degrees of freedom at bound state, (2) ‘disappearing the locked ligand’, turning off the interaction between ligand and its surroundings, (3) ‘translocating ligand’ of which corresponding free energy is zero, (4) ‘reappearance of the locked ligand’, turning on the interaction between ligand and its surroundings at bulk state, and (5) ‘unlocking the ligand’, releasing of the three restraints from [1]. Lenselink et al. tested binding free energies of congeneric ligands to four different using FEP+ from Schrödinger and obtained results in great agreement with our experimental results [27].

5. Laboratory validation

In order to validate the binding of the Drug Bank compounds, colorimetric kinetic assays were conducted. To perform these tests recombinant FBPase protein was over-expressed and purified. The construct was transformed into competent cells, and the purified construct used to over-express protein, which was harvested and purified prior to use.

5.1 Transformation(s) and purification of FBPase plasmid

Plasmids containing FBPase sequence underwent transformation with XL Blue super-competent cells [13]. Cells were plated on Luria Bertani (LB) agar plates and colonies selected for 5 mL overnights for plasmid purification using various kits [28]. Following purification, FBPase plasmids were screened for integrity and run on an agarose gel and sequenced. A transformation protocol for over-expression was then performed as previously described [13]. The contents of the tubes were transferred to LB ampicillin plates using a sterile technique. Plates were incubated for 18 h at 37°C and stored for 2.5 weeks at 4°C.

5.2 Protein over-expression, isolation, and purification

To prepare the recombinant human or closely related pig kidney FBPase enzymes for kinetic and binding assays, the recombinant proteins were over-expressed, isolated, and purified as previously described using ampicillin resistance for selection of cells containing the construct [13]. Briefly, the host cell translation was inhibited with 34 mg/mL chloramphenicol (in isopropanol). The solution was shaken again for 2–3 h at 37°C to ensure optimal growth of the host cells. Cells were isolated by pelleting in 250 mL flasks in a centrifuge at 4000 rpm. After the cells were frozen, each cell pellet was resuspended in 20 mL of 50 mM Tris pH 7.5. The supernatant was lysed via sonication to release cell contents. Each protein solution was sonicated as previously described. Sonication settings were at 10% duty cycle for 5 min, pulsing 10 s on, and 10 s off ×3. Each supernatant cell lysate was centrifuged for 30 min at 13,500 rpm at 4°C and transferred to dialysis tubing for dialysis in 50 mM Tris pH 8.0. Protein was then purified via NTA nickel affinity column as previously described [13]. In addition, gel filtration was run on a G250 column in Tris buffer pH 7.5 in 0.150 M NaCl as eluent buffer to investigate the oligomeric status of inhibitor-bound protein compared to FBPase enzyme alone.

5.3 Characterization of purified recombinant enzymes and preparation for enzyme kinetic assays

The purity of the recombinant FBPase was assessed via SDS-PAGE, and the oligomeric state was identified with native gel electrophoresis. An SDS-PAGE gel electrophoresis was used to separate and identify proteins with the correct molecular weight [9]. For the kinetic assay, a standard curve was obtained from an ammonium molybdate malachite green inorganic phosphate assay (OD660 nm) [5, 6]. The purified FBPase protein was dialyzed ×3 in 50 mM Tris buffer at pH 7.5 at 4°C. For enzyme concentration, absorbency values were recorded using spectrophotometer readings at OD280 nm. To quantify enzyme concentration, the enzyme concentration was calculated based on the quantity of micrograms per microliter (μg/μL) present. Final concentration FBPase enzyme was between 2.5 and 5.0 μg/μL. Using a colorimetric assay as described below, specific activity (SA) was determined in the absence and presence of Drug Bank inhibitors.

5.4 Kinetic assays on FBPase/ligand complexes

A colorimetric malachite green kinetic assay was utilized for FBpase activity levels to validate inhibition of the enzyme predicted with selected Drug Bank molecules. FBPase cleaves fructose 1,6-bisphosphate to fructose 6 phosphate and inorganic phosphate. The malachite green colorimetric assay is based on the change in color from brown to blue, observed when a complex is formed between malachite green, ammonium molybdate, and the product inorganic phosphate. For calculating Ki values, data was collected with varying substrate FBP concentrations (100–500 μM at pH 7.5). Malachite green dye was prepared under acidic conditions to activate the color change and quench the FBPase activity at fixed time points [13]. Absorbance readings were recorded at OD660 nm. Data collected was input to the inorganic phosphate standard curve equation to calculate product formed. Micromoles of product inorganic phosphate were calculated to determine the specific activity (SA) of FBPase enzyme +/− inhibitor. The IC50 assay was designed with fixed high concentration of substrate varying the inhibitor concentration to calculate IC50 values. For each Drug Bank ligand, the kinetic assay was performed in duplicates of triplicates. All kinetic data were fit using origin software based on methods previously described [13]. More specifically, for kinetic parameters established to fin Ki’s via curve fitting the equation y = d + (a – d) / 1 + (x / c)b was used where “x” is represented by apparent Ki when y = Vmax. IC50s were also determined with origin curve-fitting software [29].

To avoid high background due to malachite green dye interaction with some Drug Bank ligands cleavage assays utilizing phosphoglucose isomerase and glucose-6-phosphate dehydrogenase were used as coupling enzymes in validation assays for FBPase [18]. For specific activity measurements, reduction of NADP to NADPH was monitored by absorbance at 340 nm. Other assays used the same coupling enzymes but monitored the formation of NADPH by its fluorescence emission at 470 nm using an excitation wavelength of 340 nm. Assays were performed at 22°C in 50 mM Hepes, pH 7.5. Data for inhibitor (inhibition) were fit to several models using origin or sigma plot software with a Hill equation model [6]. These assays are robust and unbiased as they follow the ASBMB standard for rigor and reproducibility [30]. In addition, these assays have been published in many peer reviewed research articles [1, 2, 3, 4, 5, 14].

6. Case study analysis and discussion

6.1 Case study: allopurinol as a scaffolding molecule for drug design

The approved Drug Bank compound allopurinol is already on the market to relieve symptoms of gout. Also known as Zyloprim and Aloprim, allopurinol is composed of a dihydroxy substituted pyrimidine ring that is fused to a pyrazole ring. Allopurinol is shown below in green in Figure 3A. Allopurinol alone is not large enough to bind to the entire allosteric binding site, and more than one allopurinol molecule may bind at a time resulting in a 1:2 binding ratio. For every target allosteric binding site, 2 allopurinol compounds can fit in the allosteric binding pocket. In this study, allopurinol molecules docked in the FBPase AMP binding site were used as scaffolding molecules upon which functional groups were added to enhance their binding affinity for the FBPase. As shown below in Figure 3A, allopurinol top scoring position (Ki ~300 micomolar) overlays with the crystallographic AMP molecule adenine ring. Interestingly, with another allopurinol docked, the predicted Ki dropped to ~150 micromolar. The residues shown here Phe184, Leu173, Glu20, Val160, Arg140, Tyr113, Lys112, Glu 29, and Thr27 were all separately made flexible during reiterative runs of docking followed by full-scale MD runs for scoring predicted Ki’s. In Figure 3B, the Zinc38643891 (pink) literally straddles the AMP binding site. The predicted Ki is ~50 nM for this catechin derivative discovered in the PubChem database after similarity searches for catechin EGCG as previously described (ref). Note when compared to the crystallographic position of the AMP co-crystallized with FBPase in PDB 1FTA, the predicted position of the phenyl ring aligns in the position where the adenine ring is located. Leu174, Phe184, and Val60 alternate to interact with the zinc molecule and stabilize its position during simulations. During simulations the helix behind these three residues shifts slightly out of position causing a “domino effect” and the zinc molecule moves deeper into the AMP binding site cleft.

Figure 3.

Docking of allopurinol and Zinc38643891 into the AMP allosteric binding site of FBPase.

Figure 3C shows the best docking pose predicted for Zinc38643891 and allopurinol. This figure illustrates how allopurinol may be modified to bind tighter by derivatizing the heterocyclic ring. Details of the synthesis are preserved for a future manuscript.

Figure 3A shows the 1FTA protein coordinates with a close-up of the AMP binding site with the crystallographic ligand AMP in dark gray. The best score for the docked allopurinol is in green. Figure 3B shows the same coordinates docked with a catechin derivative (pink) from the PubChem database cross-referenced to the zinc database as Zinc38643891. Figure 3C shows the best docking results of the two compounds Zinc38643891 (pink) and allopurinol (green).

Figure 4A shows allopurinol derivative (orange) and Zinc38643891 (blue) in an FBPase activity assay showing relative activity as a function of compound concentration. The average specific activity of each concentration (triplicates) was normalized against FBPase with no inhibitor. Figure 4B shows the relative viability of HepG2 cells as a function of inhibitor concentration. The raw data from the MTT assay was normalized with HepG2 cells with no inhibitor present.

Figure 4.

Inhibition and toxicity assays comparing Zinc38643891 and sllopurinol derivative.

As shown in Figure 4A, both the allopurinol derivative and Zinc inhibit the activity of the isolated recombinant FBPase activity with ~50% inhibition in the low nanomolar range. The allopurinol derivative reaches the 50% inhibition mark at exactly 100 nM whereas the Zinc compound reaches 50% inhibition at ~25 nM, nearly 4-fold more potent than the allopurinol derivative in this cell-free activity assay. This is considered a promising hit for moving forward in the drug development pipeline to a lead compound.

However, it was surprising in Figure 4B that within the same concentration range, the Zinc38643891 was at 50% viability in the MTT viability assay with the HepG2 cells. Whereas the allopurinol derivative remained 90% viable in this range and beyond. In fact, allopurinol did not show significant drop in viability until the mid-high micromolar range. Initially, the expectation was that a derivative of a natural product catechin, the Zinc compound would have a better toxicological profile. The allopurinol derivative of Zinc38643891 was able to overcome the toxicological barrier of the original zinc molecule by substituting out the fluorophenyl rings for allopurinol in a novel synthesis protocol (manuscript under preparation). Insights from this case study that have led to the synthesis of the novel allopurinol derivative were discovered by visual analysis of the output frames of the MD runs of the allopurinol and Zinc38643891 protein-ligand complexes. During this visual inspection process (which was laborious), it became apparent that the allopurinol molecule(s) in the allosteric binding site were localizing in the same area as Zinc38643891 fluorophenyl rings. The VDDP was an integral part of the drug discovery process for this project. Currently, the lead compound(s) from this case study are being evaluated for in vivo studies.

7. Conclusions

Using our virtual drug discovery platform, symmetric computing has identified potential therapeutic small molecules used as scaffolding molecules that are approved drugs. The repurposed compounds were selected based on theoretical binding score (Ki) of the FBPase protein-compound complex. Initial validation in a colorimetric enzymatic inhibitory assay and a toxicity assay led the way to derivatives of these repurposed drugs from validated hits to lead optimization. The next steps to validate the predicted activity of each of these compounds would be an animal model (in vivo studies targeting the liver) utilizing a rat animal model. In the future, a decision will be made on which Drug Bank derivative will be advanced to human clinical studies.

Acknowledgments

We thank Jason Gao of Northeastern University for the use of the NMR facility to track synthesis of the Drug Bank derivatives. Also thank you to Jeffrey Turner and the symmetric computing support staff for this project, and Professor Nurit Haspel for valuable insight and discussions on the molecular dynamics protocols. As well as the Biotechnology Division at Roxbury Community College. Thank you to Northeastern Professional Studies College Department of Biotechnology for assistance with HepG2 toxicity assays.

Conflict of interest

None.