Novel glucokinase activators: A structure-based pharmacophore modeling, QSAR analysis, and molecular dynamics approach

Mansour Al-Sayed Ahmad; Belal O. Al-Najjar; Ashok Shakya

doi:10.3897/pharmacia.71.e131072

Abstract

Glucokinase (GK) activators are promising candidates for type 2 diabetes treatment. This study utilized structure-based pharmacophore modeling and QSAR analysis to identify novel activators. Virtual screening of a 250,000-compound library yielded eight new candidates with significant in vitro activity (over 50% activation at 25 µg/mL) and diverse structures. Molecular dynamics simulations revealed a potential mechanism involving a transient loop flip in the GK allosteric site, aligning with known activator behavior. The leading candidate, NSC12516, displayed superior hydrogen bonding with key residues (compared to known activator MRK501) and potentially stronger binding affinity due to favorable energetics. These findings clear the way for developing potent and selective GK activators for type 2 diabetes.

Keywords

glucokinase, virtual screening, pharmacophore, QSAR, molecular dynamics

Introduction

Glucokinase (GK), also known as hexokinase 4, belongs to the hexokinase family and is an inducible enzyme that is made up of 465 amino acids with a molecular mass of 52 kDa. GK’s three-dimensional structure (Fig. 1) is composed of large, small, and connected domains, with the connected domain consisting of three linker segments and serving as the primary active region. The glucose and GK binding sites are located in this area. (Matschinsky 2002; Ren et al. 2022; Su et al. 2023).

Figure 1.

Solid ribbon representation of glucokinase protein structure (PDB ID: 3F9M). Generated by Biovia Discovery Studio® Software.

GK exhibits three different forms, namely closed, open, and super-open conformations, depending on the binding of endogenous ligands or substrates to the connected domain (Fig. 2). The closed and open forms represent the states of open receiving and closed processing, respectively, and play a crucial role in determining the conversion of glucose to glucose-6-phosphate. Conversely, in the super-open form, GK is inactive and does not engage with the substrate, specifically glucose. Consequently, during instances of low blood glucose levels, the super-open form prevails, affecting the activity of GK., which is present mainly in the super-open form. (Choi et al. 2013; Toulis et al. 2020).

Figure 2.

Schematic representation of the three conformational transformations of glucokinase.

Since the 1960s, numerous studies have been conducted that have highlighted the dual functions of GK. The enzyme was identified as a glucose sensor in pancreatic β-cells and as a pacemaker in hepatic glucose conversion to glycogen. Further research discovered that GK expression in pancreatic β-cells and hepatocytes was differentially regulated by glucose and insulin, respectively. (Matschinsky 2009) (Nauck et al. 2021).

3D pharmacophores provide a highly abstract representation of the binding modes between ligands and their macromolecular targets. This abstraction allows for rationalization of binding modes even for chemically diverse ligands. Furthermore, it enables efficient virtual screening of molecular databases, facilitating rapid identification of potential ligands with desired binding characteristics. The use of 3D pharmacophores streamlines the process of identifying and evaluating new molecules with therapeutic potential. (Schaller et al. 2020) (Al-Najjar et al. 2023) For many years, researchers have utilized quantitative structure-activity relationships (QSAR) as a valuable tool in establishing connections between the physicochemical characteristics of chemical substances and their biological effects. This has enabled the development of dependable statistical models for forecasting the activities of novel chemical compounds. The fundamental concept driving this methodology is that distinctions in structural attributes play a pivotal role in accounting for the variations observed in the biological activities of these compounds.

Within this context, 3D-QSAR has emerged as a logical progression from traditional methods. It leverages the three-dimensional attributes of ligands to anticipate their biological effects by employing robust chemometric techniques. (Al-Anazi et al. 2022) (Verma et al. 2010).

Since the initial introduction of the first glucokinase activator in 2003, numerous activators belonging to this class have been developed and subjected to evaluation. These activators are characterized by their small molecular size, enabling them to bind to a specific allosteric site on the enzyme. When these polymers occupy this unique site, they promote the stabilization of a conformation of GK that exhibits a notably high affinity, ultimately facilitating the activation of the enzyme. Remarkably, it has been observed that this allosteric cleft is in the hinge region (Fig. 3). (Sarabu et al. 2011).

Figure 3.

A solid ribbon representation of the glucokinase protein showing the location of the allosteric site in the GK enzyme. Generated by Biovia Discovery Studio®.

Many chemical groups were reported by previous studies as potential GK activators, like azole derivatives, pyridine derivatives, azole-pyridine derivatives, or miscellaneous. (Grimsby et al. 2003; McKerrecher et al. 2006; Mitsuya et al. 2009).

But the use of older-generation GKAs raised concerns about both efficacy and safety. Key issues included the risk of hypoglycemia, induction of fatty liver disease, dyslipidemia, and a decrease in long-term efficacy. The appearance of hypoglycemia and dyslipidemia, which were attributed to excessive stimulation of pancreatic and liver glucokinase, respectively, were identified as potential risks early in the development of GKA (Matschinsky 2013; Agius 2014). That is why we continue researching for new, safe, and effective GKAs.

Experimental

In-silico study

In this study, we adopted the structure-based pharmacophore modeling procedure reported by Al-Najjar and coworkers (2023). The study initially started by generating valid pharmacophore models using Biovia Discovery Studio 3.1 (Biovia Dassault Systems, Discovery Studio 2016). Those models were subjected to QSAR analysis to generate a representative equation that explains the activity. Next, the selected pharmacophore models that appeared in the QSAR equation were used to screen 264,000 compounds from the NCI database. Finally, the successful candidates were tested for their biological activity. (Al-Najjar et al. 2023).

Selecting the crystal structure

More than 50 X-ray crystal structures in the Protein Data Bank (https://www.rcsb.org) for glucokinase were reviewed to choose one structure according to the following criteria: (1) high resolution; (2) without glucokinase regulatory protein in its structure; and (3) co-crystallized with an activator.

The crystal structure of 3F9M co-crystallized with its activator, MRK501, was selected according to the mentioned criteria.

Validation data set

The generated pharmacophore models were validated by assessing their abilities to selectively capture diverse active compounds from a large testing list of actives and decoys. The list of studied compounds as glucokinase activators was downloaded from the ChEMBL database (https://www.ebi.ac.uk/chembl/).

The list was composed of synthetic or organic compounds listed as glucokinase activators and then divided into active and inactive compounds according to their EC50 values to make a training set that serves as a template for evaluating the predictive accuracy of pharmacophore models in identifying known active compounds.

Pharmacophore model generation

The data set prepared previously was uploaded to the ‘Pharmacophore Generation Protocol’ impeded in Biovia Discovery Studio® software. The inputs in the protocol were the (3F9M) protein and the training set, and then we changed a number of parameters to create the different pharmacophoric hypotheses. The parameters were (omitted features between 0 and 1; minimum interference distance 2, 2.5, and 3 Å; and maximum exclusion volume 4 and 5 Å3); and the output is a set of pharmacophores suggested to be suitable for interaction with the activators at the allosteric site. An internal validation was performed by supplying a collection of both active and inactive ligands, and the outcomes were visualized in the form of a receiver operating characteristic (ROC) curve.

Pharmacophore mapping

More than 100 pharmacophores were generated in step (3), and then the pharmacophores were tested against a test data set to choose the best pharmacophore that could recognize the active compounds and differentiate them from the inactive compounds.

QSAR analysis

The selected pharmacophores were employed in a ligand pharmacophore fitting approach for the ligands, using the optimal fit option in DS. The calculated descriptors contained (log P, molecular weight, number of aromatic numberings, number of H-bond acceptors, number of H-bond doners, number of rotatable bonds, polar surface area, extended-connectivity fingerprints (ECFPs), as well as the fitting values of the generated pharmacophores). The results were uploaded to the ‘Create Genetic Function Approximation Model’ protocol.

Virtual screening

The QSAR equation and the chosen pharmacophore models were used as parameters to virtually screen a large chemical database composed of more than 250,000 compounds obtained from the National Cancer Institute (https://www.cancer.gov/).

A total of 250 conformers were generated using the BEST conformation method, ensuring comprehensive and enhanced coverage of the conformational space. To be considered a hit, a selected molecule from the screening process should successfully match all the pharmacophoric features. The fit values and associated molecular descriptors of each hit were incorporated into the previously prepared QSAR equation. The top-ranked compounds, determined through QSAR predictions, were subsequently subjected to in vitro evaluation against the glucokinase enzyme.

In vitro assay

Taha et al. (2015) reported an in vitro assay procedure that was adopted for this study.

The bioassay is based on the enzymatic phosphorylation of D-glucose by GK, which results in the formation of D-glucose-6-phosphate. This product is further oxidized by the glucose-6-phosphate dehydrogenase (G6PD) enzyme in the presence of NADP, leading to the production of 6-phospho-D-gluconate and NADPH. The NADPH molecule exhibits maximum absorbance (λmax) at 340 nm. The rate at which NADPH is generated is directly proportional to the catalytic activity of GK. (Taha et al. 2015).

In the experimental procedure, stock solutions of the test samples are initially prepared using dimethyl sulfoxide (DMSO). These stock solutions are then diluted in deionized water in a serial manner to achieve the desired working concentrations for the bioassay.

All chemicals needed for the bioassay were acquired from Sigma-Aldrich Company and used without further purification.

To perform the bioassay, the following procedure was followed:

A total of 10 µL of the tested sample solution was added to a reaction mixture containing the following components: Tris-HCl buffer (75 mM, pH 9.0 at 30 °C) -24 mL, MgCl2 (600 mM in deionized water) -1 mL (equivalent to 20.10 mM in the reaction mixture), ATP (120 mM in deionized water) -1 mL (equivalent to 4.02 mM in the reaction mixture), D(+)glucose (360 mM in deionized water) -1 mL (equivalent to 12.10 mM in the reaction mixture, which is close to the S0.5 reported for GK, i.e., 8 mM), NADP (27 mM in deionized water) -1 mL (equivalent to 0.90 mM in the reaction mixture). Subsequently, G6PD (1000 U/mL in cold deionized water) -3 µL (equivalent to 0.031 U/L in the reaction mixture) was added. This was followed by the addition of human GK solution (Sigma-Aldrich, CAS Number 9001-51-8) (0.05 units/mL) -3 µL (equivalent to 1.56 × 10^-6 U/L in the reaction mixture), which was prepared in cold tris buffer (pH 8.5, 4 °C). The addition of the human GK solution initiated the reaction in the bioassay (in all experiments and controls, samples were prepared in triplicate).

Molecular dynamics simulation

Molecular dynamics simulations were initiated and conducted for 25 ns.

Minimization: This work employed the PMEMD module to optimize the structure of protein-ligand interactions. There were three steps to energy minimization. The first step involved relaxing the water molecules surrounding the complex over 5,000 steepest descent steps, followed by another 5,000 conjugate gradient steps. The second phase sought to eliminate any collisions between the side chains for the same number of cycles. Finally, the entire system was minimized for 10,000 conjugate gradient cycles more. During minimization, the non-bonded contact cut-off was set at 9.0 Å, with constant volume periodic borders and shake applied exclusively for hydrogen atoms.
Heating: The following phase involved heating the systems progressively over a period of 6 ns in three stages to begin running the molecular dynamics simulations with the assistance of the PMEMD module. Heating in many steps allows the system to equilibrate at each desired temperature, preventing the simulated system from breaking down. In this study, the systems’ temperature was raised from 0 to 100 °K in the first phase, 100 to 200 °K in the second step, and 200 to 310 °K in the third step utilizing an NVT ensemble.
Equilibration: The systems must be balanced in order to maintain the temperature (NPT ensemble) and pressure (isotropic position scaling) constant. Three stages of six ns each were employed to equilibrate the systems using PMEMD. Throughout the equilibration, 3 and 1 kcal (mol. Å2)-1 constraint forces were used for the first and second steps, respectively. A cut-off for non-bonded interactions was established at 9.0 Å.
Production: Running a production simulation for 25 ns at 310 °K was the last step in the molecular dynamics’ simulations. No constraint force was used when using the NPT ensemble. 9.0 Å was chosen as the cut-off for non-bonded interactions. Root mean square deviation (RMSD), root mean square fluctuation (RMSF), and hydrogen bond interactions were examined on the trajectories using MDAnalysis and CPPTRAJ. (Al-Najjar and Saqallah 2023) (Case et al. 2023).

Results and discussion of pharmacophore models and validation

The receptor-ligand pharmacophore generation protocol in Discovery Studio was used to create more than 100 pharmacophore models. The models were generated by varying parameters such as maximum features, minimum inter-feature distance, and maximum exclusion volume distance.

To assess the efficacy of the pharmacophoric model in distinguishing active and inactive compounds, the area under the ROC curve was calculated. (Suppl. material 1) presents the results of the ROC analysis for the selected pharmacophores, which will be used in the subsequent QSAR analysis. The pharmacophore models primarily consist of hydrogen bond donor (HBD), aromatic ring (RingArom), and hydrophobic (Hbic) features. Some models also incorporate additional excluded volumes to account for spatial constraints.

QSAR analysis

In this study, classical QSAR analysis was performed to identify the optimal combination of pharmacophores and other 2D and 3D descriptors that can explain the biological activities of the ligands. The goal was to develop a robust QSAR equation that accurately predicts the activities of novel compounds.

Multiple linear regression and GFA (genetic function approximation) were employed to obtain the best QSAR equation. The fit values obtained from the pharmacophoric mapping of the dataset were combined with various physicochemical descriptors of the ligands.

By integrating the information from the pharmacophoric models and the additional descriptors, the QSAR analysis aimed to establish a quantitative relationship between the structural and physicochemical properties of the ligands and their biological activities. This approach allows for the identification of key molecular features and properties that contribute to the observed activity, facilitating the design and discovery of novel pharmacologically active compounds. The selected QSAR equation is:

Log (1/EC50) =-1.0397 + 0.56233 * ALogP + 3.5554 * Count<ECFP_6:914325265>+ 1.3589 * Count<ECFP_6:-797085356> - 1.1095 * Count<ECFP_6:-177264675>+ 0.75328 * Count<ECFP_6:-1059365320> + 3.7999 * Count<ECFP_6:864909220> - 1.997* Count<ECFP_6:432225086> + 0.53069 * R1T1D2.5P6 - 0.23105 * R2T1D3P10 + 0.43508 * R2T1D3P7

r ² = 0.9864 r²(adj) = 0.9767 r²(pred) = 0.9180

Equation 1; where EC50: Half maximal effective concentration; ALogP: Log of the octanol-water partition coefficient using Ghose and Crippen’s method; ECFP: Extended-Connectivity Fingerprints; r2 is the correlation coefficient against the testing compounds; r2(adj) is r2 adjusted for the number of terms in the model; r2(pred) is the prediction r2, equivalent to q2 from a leave-1-out cross-validation.

In the QSAR equation (Equation (1)), R1T1D2.5P6 refers to the fit values of the ligands against the sixth pharmacophoric hypothesis. This calculation involves the interfacial distance of 2.5 Å from the first trial and so on for the other two pharmacophoric hypotheses present in the equation (R2T1D3P7 and R2T1D3P10). The inclusion of the pharmacophoric model in the equation indicates that two of the pharmacophores (which are the pharmacophoric hypotheses with the positive factors R1T1D2.5P6 and R2T1D3P7) are able to effectively explain the observed bioactivities. But on the contrary, if the ligands fit the third pharmacophoric hypothesis (R2T1D3P10), it will be bad for activity.

The QSAR equation and its associated parameters are utilized to establish a quantitative relationship between the structural features of molecules and their biological activities, so it becomes possible to predict and analyze the bioactivity of similar ligands based on their adherence to the defined pharmacophoric model.

Fig. 4 shows the pharmacophoric model, which is composed of HBD, RingArom, and Hbic features mapped to the co-crystallized ligand.

Figure 4.

A. R1T1D2.5P6 pharmacophore model; B. 3D representations of the chemical structure of the co-crystallized ligand (MRK501) obtained from crystal structure 3F9M mapped to R1T1D2.5P6 phramacophore, wherein orange: aromatic ring, cyan: hydrophobic, and magenta: hydrogen bond donor features.

Virtual screening

As mentioned earlier, a data set consisting of 250,000 compounds obtained from NCI was screened against the QSAR equation, and the pharmacophore models with the highest ranking of 30 compounds based on their predicted EC50 (Suppl. material 2) were selected for subsequent in vitro analysis.

In vitro assay

During the in-vitro assay of selected compounds, absorbance readings were recorded. Out of the tested compounds, eight exhibited significantly higher positive activity compared to the control (the reaction mixture without any tested compound). The chemical structures of the specific compounds demonstrating such activity are detailed in Suppl. material 3.

Each one of the eight compounds was tested at (25 μg/mL) concentration, and the test was repeated three times. Then the average reading of each compound was divided by the average reading of the vehicle. The results are summarized in Fig. 5.

Figure 5.

Graphical summary of in-vitro results.

We noticed that at 25 µg/mL concentration, all the compounds showed good activity except (10666), which demonstrated moderate activity as a GK activator.

As we can see, the chemical structures for the eight compounds are miscellaneous with no similarity to any reported group, but in silico they preserve the interactions with the pharmacophore.

Suppl. material 4 shows the chemical structures of the tested compounds and compares their binding in the protein crystal structure (3F9M) with how it maps against the R1T1D2.5P6 model without conformational adjustments. The main intermolecular interactions performed by the co-crystal structure MRK501 (Suppl. material 4: fig. S1A) were hydrogen bond interactions with Tyr215 and Pi-anion interactions with Glu221 and hydrophobic interactions with Val62, Pro66, Ile159, Ile211, Met235, and Pi-Pi interactions with Tyr214. Similar amino acids were involved in the interactions with the 8 NCI compounds, as seen in Suppl. material 4: fig. S1B–I. Similar amino acids that were found to interact with the novel compound, implies that the compound may adopt a similar binding mode and exploit the same or similar binding interactions as the known active compound (MRK501), which indicates the potential for similar activity and specificity.

In QSAR analysis, parameters appeared in the equation, i.e., ALogP: Log of the octanol-water partition coefficient; ECFP: extended-connectivity fingerprints; and the pharmacophoric model R1T1D2.5P6. This indicates a positive correlation between those parameters and the activity predicted. The positive correlation of AlogP may indicate that compounds with higher LogP values tend to have a greater ability to interact with hydrophobic binding pockets, resulting in increased binding affinity and activity. The positive correlation between a specific ECFP and the predicted activity means that the presence or occurrence of that particular ECFP feature tends to be associated with higher activity levels. In other words, compounds that possess this ECFP feature are more likely to exhibit the desired activity.

Finally, experimental validation of the previous in silico procedure has shown that selected compounds were active, which provides strong evidence that the predicted pharmacophore models and the pharmacophoric features are functionally relevant. The activity of the selected compounds indicates that these interactions play a crucial role in the binding and potential efficacy of the compound.

Molecular dynamics analysis

To validate our in vitro bioassay findings and gain deeper insights into the compound-glucokinase (GK) interactions, we performed 25 ns molecular dynamics (MD) simulations. (Hariono et al. 2019)

The first system simulated the crystal structure of GK complexed with the known activator MRK501 (PDB ID: 3F9M). The second system investigated NSC12516, the compound exhibiting the highest in vitro activity. The root mean square deviation (RMSD) of the protein backbone over the simulation trajectory was monitored (Fig. 6A). Both systems demonstrated relative stability, with average RMSD values of 1.66 ± 0.23 Å for NSC12516 and 1.69 ± 0.21 Å for MRK501, indicating minimal structural drift during the simulations.

Figure 6.

Visual representation of the Molecular Dynamics Simulation Analysis, where A. is root square mean deviation (RMSD) and B. is root mean square fluctuation (RMSF); blue: NSC12516; red: MRK501.

Detailed analysis of the RMSF profile (Fig. 6B) revealed a transient increase in fluctuations between 15 and 16 ns. This corresponds to a loop flip involving Pro66, located in the allosteric site known to interact with activators. Interestingly, this observation aligns with previous findings on GK activator-induced conformational changes. (Liu et al. 2020).

These results suggest that the identified activators might induce a similar loop movement, potentially contributing to their mechanism of action.

The observed fluctuations around residues 50–65 suggest potential movement within the broader allosteric region beyond the immediate binding pocket. These fluctuations might be independent of ligand binding or could represent a “domino effect” propagating through the protein structure. Next, we investigated hydrogen bond interactions between the amino acid residues lining the allosteric site of GK and the bound ligands (Table 1). Fig. 7 summarizes the overall number of hydrogen bonds formed in each system. These analyses aimed to identify specific interactions potentially contributing to ligand binding and activation.

Figure 7.

Two-dimensional plots of the number of hydrogen bond interactions in the simulated systems, where A. is NSC12516 and B. is MRK501.

Table 1.

Download as

CSV

XLSX

Hydrogen bond interaction analyses of NSC12516 and MRK501 in complex with GK.

	Acceptor	Doner	Occupancies (%)	Average Distance (Å)	Average Angle
NSC12516-GK Complex		NSC12516@H43/N5	74.57	2.85	157.91
		NSC12516@H48/O2	18.61	2.71	155.38
	GLU67@OE1	NSC12516@H48/O2	17.6	2.71	155.58
	NSC12516@O2	GLU67@H/N	1.54	2.94	155.21
MRK501-Gk Complex	MRK501@N4	ARG63@H/N	4.22	2.94	154.39
MRK501-Gk Complex	MRK501@N17	SER64@HG/OG	3.28	2.91	151.58

NSC12516 was able to maintain 3 hydrogen bonds as an accepter from Arg63 (74.57%), Glu67 (81.61% and 17.6%), and 1 hydrogen bond as a doner for Glu67 (1.54%), while MRK501 maintained only 2 hydrogen bonds as an accepter from Arg63 (4.22%), and Ser64 (3.28%). Overall, the higher number and occupancy rates of hydrogen bonds formed by NSC12516 suggest that it might have a stronger affinity and potentially better stability when bound to glucokinase compared to MRK501.

To further understand the binding interactions, we calculated the energetic components between NSC12516 and MRK501 in complex with GK using the MM-PBSA method (Table 2). This analysis employed 50 frames extracted from the final 5 ns of each simulation. Interestingly, NSC12516 exhibited binding predominantly driven by non-polar interactions. Notably, the key parameter, the difference in binding free energy (ΔGbind) between NSC12516 and MRK501, revealed a lower value for NSC12516. This lower ΔGbind suggests a potentially higher binding affinity for NSC12516 compared to MRK501.

Table 2.

Download as

CSV

XLSX

Predicted energy components and binding affinities of NSC12516 and MRK501 towards the allosteric site of GK using the MM-PBSA approach for the last 5 ns of the molecular dynamics’s simulations.

Energy component (kcal/mol)	Energy value (kcal/mol) per complex
Energy component (kcal/mol)	NSC12516-GK	MRK501-GK
van der Waals Energy (∆E_vdW)	-61.9288	-25.5885
Electrostatic Energy (∆E_EL)	-34.4406	-44.1295
Polar Solvation Energy (∆E_PB)	48.1815	49.8056
Non-Polar Solvation Energy (∆E_NPOLAR)	-4.036	-2.5816
Total Gas Phase Free Energy (∆G_gas)	-96.3694	-69.7181
Total Solvation Free Energy (∆G_solv)	44.1455	47.224
Total Energy (∆G_bind)	-52.2239	-22.4941

Conclusions

This study successfully identified eight novel glucokinase activators with promising in vitro activity through a structure-based approach. Molecular dynamics simulations suggest a potential mechanism involving a loop flip in the allosteric site, and the leading candidate, NSC12516, exhibits superior binding characteristics compared to a known activator. These findings lay the groundwork for further development of potent and selective GK activators as a therapeutic strategy for type 2 diabetes. Further investigation is warranted to fully elucidate the mechanism of action and translate these promising candidates into clinical applications.

Acknowledgments

The authors would like to acknowledge Al Ahliyya Amman University for funding this project.

References

Agius L (2014) Lessons from glucokinase activators: the problem of declining efficacy. Expert Opinion on Therapeutic Patents 24: 1155–1159. https://doi.org/10.1517/13543776.2014.965680

Al-Anazi M, Khairuddean M, Al-Najjar BO, Alidmat MM, Kamal NNSNM, Muhamad M (2022) Synthesis, anticancer activity and docking studies of pyrazoline and pyrimidine derivatives as potential epidermal growth factor receptor (EGFR) inhibitors. Arabian Journal of Chemistry 15: 103864. https://doi.org/10.1016/j.arabjc.2022.103864

Al-Najjar BO, Abbas MA, Sibai OA, Saqallah FG, Al-Kabariti AY (2023) QSAR, structure- based pharmacophore modelling and biological evaluation of novel platelet ADP receptor (P2Y 12) antagonist. RSC Medicinal Chemistry 14: 239–246. https://doi.org/10.1039/D2MD00285J

Al-Najjar BO, Saqallah FG (2023) Unlocking the Conformational Changes of P2Y12: Exploring an Acridinone Compound’s Effect on Receptor Activity and Conformation. Molecules 28: 3878. https://doi.org/10.3390/molecules28093878

Case DA, Aktulga HM, Belfon K, Cerutti DS, Cisneros GA, Cruzeiro VWD, Forouzesh N, Giese TJ, Götz AW, Gohlke H (2023) AmberTools. Journal of chemical information and modeling 63: 6183–6191. https://doi.org/10.1021/acs.jcim.3c01153

Choi JM, Seo M-H, Kyeong H-H, Kim E, Kim H-S (2013) Molecular basis for the role of glucokinase regulatory protein as the allosteric switch for glucokinase. Proceedings of the National Academy of Sciences 110: 10171–10176. https://doi.org/10.1073/pnas.1300457110

Grimsby J, Sarabu R, Corbett WL, Haynes N-E, Bizzarro FT, Coffey JW, Guertin KR, Hilliard DW, Kester RF, Mahaney PE (2003) Allosteric activators of glucokinase: potential role in diabetes therapy. Science 301: 370–373. https://doi.org/10.1126/science.1084073

Hariono M, Nuwarda RF, Yusuf M, Rollando R, Jenie RI, Al-Najjar B, Julianus J, Putra KC, Nugroho ES, Wisnumurti YK (2019) Arylamide as potential selective inhibitor for matrix metalloproteinase 9 (MMP9): design, synthesis, biological evaluation, and molecular modeling. Journal of Chemical Information and Modeling 60: 349–359. https://doi.org/10.1021/acs.jcim.9b00630

Liu W, Yao C, Shang Q, Liu Y, Liu C, Meng F (2020) Insights into the binding of dorzagliatin with glucokinase: A molecular dynamics simulation. Journal of Theoretical and Computational Chemistry 19: 2050027. https://doi.org/10.1142/S0219633620500273

Matschinsky FM (2002) Regulation of pancreatic β-cell glucokinase: from basics to therapeutics. Diabetes 51: S394–S404. https://doi.org/10.2337/diabetes.51.2007.S394

Matschinsky FM (2009) Assessing the potential of glucokinase activators in diabetes therapy. Nature Reviews Drug Discovery 8: 399–416. https://doi.org/10.1038/nrd2850

Matschinsky FM (2013) GKAs for diabetes therapy: why no clinically useful drug after two decades of trying? Trends in Pharmacological Sciences 34: 90–99. https://doi.org/10.1016/j.tips.2012.11.007

McKerrecher D, Allen JV, Caulkett PW, Donald CS, Fenwick ML, Grange E, Johnson KM, Johnstone C, Jones CD, Pike KG (2006) Design of a potent, soluble glucokinase activator with excellent in vivo efficacy. Bioorganic & Medicinal Chemistry Letters 16: 2705–2709. https://doi.org/10.1016/j.bmcl.2006.02.022

Mitsuya M, Kamata K, Bamba M, Watanabe H, Sasaki Y, Sasaki K, Ohyama S, Hosaka H, Nagata Y, Eiki J-i (2009) Discovery of novel 3, 6-disubstituted 2-pyridinecarboxamide derivatives as GK activators. Bioorganic & Medicinal Chemistry Letters 19: 2718–2721. https://doi.org/10.1016/j.bmcl.2009.03.137

Nauck MA, Wefers J, Meier JJ (2021) Treatment of type 2 diabetes: challenges, hopes, and anticipated successes. The Lancet Diabetes & Endocrinology 9: 525–544. https://doi.org/10.1016/S2213-8587(21)00113-3

Ren Y, Li L, Wan L, Huang Y, Cao S (2022) Glucokinase as an emerging anti-diabetes target and recent progress in the development of its agonists. Journal of Enzyme Inhibition and Medicinal Chemistry 37: 606–615. https://doi.org/10.1080/14756366.2021.2025362

Sarabu R, Berthel SJ, Kester RF, Tilley JW (2011) Novel glucokinase activators: a patent review (2008–2010). Expert Opinion on Therapeutic Patents 21: 13–33. https://doi.org/10.1517/13543776.2011.542413

Schaller D, Šribar D, Noonan T, Deng L, Nguyen TN, Pach S, Machalz D, Bermudez M, Wolber G (2020) Next generation 3D pharmacophore modeling. Wiley Interdisciplinary Reviews: Computational Molecular Science 10: e1468. https://doi.org/10.1002/wcms.1468

Su J, Luo Y, Hu S, Tang L, Ouyang S (2023) Advances in research on type 2 diabetes mellitus targets and therapeutic agents. International Journal of Molecular Sciences 24: 13381. https://doi.org/10.3390/ijms241713381

Taha MO, Habash M, Ma’mon MH, Abdelazeem AH, Qandil A (2015) Ligand-based modeling followed by in vitro bioassay yielded new potent glucokinase activators. Journal of Molecular Graphics and Modelling 56: 91–102. https://doi.org/10.1016/j.jmgm.2014.12.003

Toulis KA, Nirantharakumar K, Pourzitaki C, Barnett AH, Tahrani AA (2020) Glucokinase activators for type 2 diabetes: challenges and future developments. Drugs 80: 467–475. https://doi.org/10.1007/s40265-020-01278-z

Verma J, Khedkar VM, Coutinho EC (2010) 3D-QSAR in drug design-a review. Current topics in medicinal chemistry 10: 95–115. https://doi.org/10.2174/156802610790232260

Supplementary materials

Supplementary material 1

Pharmacophore models generated by Receptor-ligand pharmacophore generation protocol embedded in Discovery Studio