ABSTRACT

Human sodium-glucose cotransporter 2 (hSGLT2) is a membrane protein responsible for glucose reabsorption from the glomerular filtrate in the proximal tubule. Inhibition of hSGLT2 has been regarded as a brand new therapeutic approach for the treatment of type 2 diabetes mellitus (T2DM) due to its non-insulin related characteristics with less side effects. Current commercially available hSGLT2 inhibitors are all C-glycoside inhibitors. Previous studies have reported that N-glycoside inhibitors have better potential to serve as new drugs due to their good metabolic stability. In addition, non-glycoside inhibitors have been shown to exhibit the capability to overcome the existing problems of current glycoside inhibitors, including low tissue permeability, poor stability and short serum half-time. Here, we aimed to discover novel N-glycoside and non-glycoside hSGLT2 inhibitors by a combination of several computational approaches. A ligand-based pharmacophore model was generated, well validated and subsequently utilized as a 3D query to identify novel hSGLT2 inhibitors from National Cancer Institute (NCI) and Traditional Chinese Medicine (TCM) databases. Finally, one N-glycoside (NSC679207) and one non-glycoside (TCM_Piperenol_A) hSGLT2 inhibitors were successfully identified, which were proven to exhibit excellent binding affinities, pharmacokinetic properties and less toxicity than the commercially available hSGLT2 inhibitor, canagliflozin, via molecular docking, ADMET prediction, molecular dynamics (MD) simulations and binding free energy calculations. All together, our results strongly suggest that these two compounds have great potential to serve as novel hSGLT2 inhibitors for the treatment of T2DM and their efficacies may be further examined by a series of in vitro and/or in vivo bioassays.

*Equal contribution.

Keywords:

Human Sodium-Glucose Cotransports 2 (hSGLT2), Type 2 Diabetes Mellitus (T2DM), Ligand-Based Pharmacophore Model, Molecular Docking, Molecular Dynamics (MD) Simulations

1. Introduction

Diabetes mellitus (DM) has become a severe health problem worldwide and its prevalence has rapidly increased over the past decades [1] . DM is a chronic metabolic disease characterized by the disorder of glucose homeostasis, which further causes long-term hyperglycemia and the development of high-risk vascular complications, including heart failure, nephropathy, retinopathy and neuropathy [2] . There are two common types of DM: type 1 diabetes mellitus (T1DM) and type 2 diabetes mellitus (T2DM). T1DM results from the absolute lack of insulin secretion due to its complete β-cell dysfunction, while T2DM is caused by insulin deficiency and insulin resistance [3] . Nowadays, T2DM accounts for almost 90% - 95% of DM patients and is extending towards a younger population [4] . So far, several commercially available therapeutic agents, such as biguanides, sulfonylureas, thiazolidinediones, DPP-IV inhibitors and GLP-1 receptor agonist, have been proven to exhibit hypoglycemic effects for treating T2DM by increasing insulin secretion or insulin sensitivity [5] . However, not all T2DM patients can be effectively treated by current therapies and these therapeutic agents are also limited by several adverse side effects, including weight gain, hypoglycemia, gastrointestinal flatulence, diarrhea, edema and lactic acidosis [6] [7] . Therefore, development of new therapeutic agents with novel hypoglycemic mechanism and fewer side effects is highly desired for the treatment of T2DM.

Human sodium-dependent glucose co-transporters (hSGLTs) are a group of membrane proteins responsible for glucose reabsorption from the glomerular filtrate in the proximal convoluted tubule. The mechanism by which hSGLTs mediate glucose reabsorption in the kidney is driven by the electrochemical sodium gradient between the glomerular filtrate and the proximal convoluted tubule [8] . Among them, hSGLT2, a low-affinity and high-capacity transporter expressed exclusively in the top third segment of the proximal convoluted tubule of the kidney, is responsible for 90% of renal glucose reabsorption [9] [10] . Inhibition of hSGLT2 induces glucose excretion from proximal convoluted tubule into urine and thereby reduces blood glucose concentration. Besides, hSGLT2 inhibitors have been confirmed to reduce body weight and prevent hypoglycemia [11] . Since hSGLT2 is mainly distributed in the kidney, hSGLT2 inhibitors would not cause common discomforts of other T2DM therapeutics agents, such as gastrointestinal flatulence, diarrhea and edema [12] . Moreover, recent studies have shown that hSGLT2 inhibitors also exhibit the capability to prevent patients from serious complications, including heart failure, nephropathy and liver disease [13] [14] [15] . Therefore, hSGLT2 inhibitors have been considered as effective therapeutic agents with less side effects, thus can be served as monotherapy or combined with other hypoglycemic agents due to their distinctive hypoglycemic mechanism.

Phlorizin, an O-glucoside natural product substance, is known to be the first hSGLT2 inhibitor [16] . Although phlorizin has sufficient efficacy against hSGLT2, it is ultimately considered to be inappropriate as a hypoglycemic agent due to gastrointestinal side effects, low oral bioavailability and insufficient selectivity for hSGLT2. Subsequently, several O-glycoside inhibitors were developed to increase the selectivity for hSGLT2 but the shortcomings were not significantly improved [17] . Due to the fact that C-glycoside inhibitors are more metabolically stable because of their resistance to gastrointestinal β-glucosidases, they exhibit better efficacy and have fewer side effects comparing to O-glycoside inhibitors [18] [19] [20] . So far, three C-glycoside inhibitors, canagliflozin, dapagliflozin and empagliflozin, have been approved by the U.S. Food and Drug Administration (FDA) for the treatment of T2DM [21] . Although C-glycoside inhibitors exhibit sufficient efficacy, previous studies indicated that N-glycoside inhibitors may exhibit better potential to serve as a drug target for the treatment of T2DM due to their good metabolic stability, oral bioavailability, and low clearance [22] [23] . Comparing to glycoside inhibitors, non-glycoside inhibitors show several advantages, such as better tissue permeability, stability, serum half-time, less expensive and easier to synthesize [24] [25] . Therefore, development of novel N-glycoside and non-glycoside hSGLT2 inhibitors would provide greater advance in antidiabetic therapy.

Over the past few years, since the crystal structure of hSGLT2 is not yet available, the strategies for developing new hSGLT2 inhibitors are limited to chemical synthesis, structure-activity relationships (SAR) analysis, and lead modification [26] [27] [28] . These strategies would be high cost and time-consuming in the process of new drug development and the results are also limited to the highly structural similarity. Thus, in this study, we aimed to discover novel N-glycoside and non-glycoside hSGLT2 inhibitors through a combination of several computational approaches, such as ligand-based pharmacophore modeling, virtual screening, homology modeling, molecular docking, ADMET predictions, molecular dynamics (MD) simulations and binding free energy calculations. These computational approaches have been successfully and intensively used for the rapid assessment of small molecular databases in order to lower the cost and speed up the early stage of the new drug development process. Finally, one N-glycoside (NSC679207) and one non-glycoside (TCM_Piperenol_A) hSGLT2 inhibitors successfully identified from National Cancer Institute (NCI) and Traditional Chinese Medicine (TCM) databases, respectively, were shown to exhibit better binding affinities, pharmacokinetic properties and less toxicity comparing to the currently commercially available canagliflozin and can be further examined by a series of in vitro and/or in vivo bioassays to ensure their efficacy. Besides, our results also reveal several essential interactions between these two compounds and the binding site of hSGLT2, which may ultimately contribute to the design of novel and potent hSGLT2 inhibitors for the clinical use in the future.

2. Materials and Methods

2.1. Preparation of the Training and Test Set Compounds

A total of 45 known hSGLT2 inhibitors were collected from the previous literatures [29] - [36] . Among them, 22 compounds, whose IC₅₀ values range from 0.3 to 4600 nM, were selected as the training set (Figure 1(A)); while the remaining 23 compounds, whose IC₅₀ values range from 2.39 to 1376 nM, were selected as the test set (Figure 1(B)). The 2D chemical structures of these compounds were sketched by ChemSketch and then converted into 3D structures by BIOVIA Discovery Studio (DS) 2017 R2. These compounds were subsequently optimized under the ionization state at pH 7.0 using the Prepare Ligand protocol of DS 2017 R2.

2.2. Ligand-Based Pharmacophore Model Generation

Before starting the pharmacophore model generation process, the conformers of each training set compound were generated by the BEST conformation model generation method and limited to a maximum conformer value of 255 with a 20 kcal/mol energy cut off. The multiple acceptable 3D conformations were obtained and then were inputted into DS 2017 R2 to generate 10 possible pharmacophore hypotheses using HypoGen module, which include hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), hydrophobic (HY), ring aromatic (RA) and excluded volumes (EV). The chemical features were set to a minimum of 1 to a maximum of 5 and the maximum of EV was set to 5, while the other parameters were set to default values [37] . These generated hypotheses were ranked according to their total cost, correlation coefficient, and root mean square deviation (RMSD). In order to discriminate the quality of these hypotheses, two theoretical cost values, null cost and fixed cost, were evaluated. The null cost represents the highest cost of the model with no correlation, while the fixed cost represents the cost of the simplest model which fits all data perfectly. According to the above mentioned criteria, the hypothesis with the lowest total cost, the lowest RMSD value, and the highest correlation coefficient was considered as the best pharmacophore model.

2.3. Ligand-Based Pharmacophore Model Validation

The quality of the best pharmacophore model was validated by test set prediction and Güner-Henry (GH) score method [38] . 23 test set compounds were used to elucidate whether the generated pharmacophore model is proficient to predict the activities of the compounds other than the training set compounds. GH score method was used to determine how well the generated pharmacophore model could differentiate potential hSGLT2 inhibitors from other non hSGLT2

inhibitors. A total of 550 compounds, including 50 active compounds selected from the previous literatures [29] - [36] and 500 decoys from DUD-E website, were used to calculate the GH score. These 50 active compounds were different from the 45 active compounds in the training and test sets. All compounds were used to construct the 3D database (maximum of 255 conformations for each compound) using Build 3D Database of DS 2017 R2. Then, the generated pharmacophore model was employed for screening the 3D database using Screen Library with default parameters. A set of statistical parameters were obtained and then the GH score was calculated to evaluate the quality of the generated pharmacophore model.

2.4. Virtual Screening

The overall flowchart showing the entire virtual screening process and the screening criteria in this study areee presented in Figure 2. First, the NCI and TCM databases were filtered by the Lipinski rule of five [39] and Veber rule [40] of DS 2017 R2 to select the compounds with good drug-like properties. Then, the well validated pharmacophore model was utilized as a 3D query to identify hSGLT2 inhibitors using Screen Library with default parameters. Compounds that successfully mapped all the features of the pharmacophore model were considered as potential hit compounds, which were further submitted to molecular docking. The hit compounds with higher binding affinities comparing to that of canagliflozin, the commercially available hSGLT2 inhibitor, were considered as novel hSGLT2 inhibitors. Finally, the pharmacokinetic properties and binding stabilities

of these novel hSGLT2 inhibitors were further evaluated through ADMET prediction and MD simulations, respectively.

2.5. Homology Modeling

The amino acid sequence of hSGLT2 (accession code: NP_003032) was obtained from the National Center for Biotechnology Information (NCBI) protein sequence data bank. Since the crystal structure of vSGLT (PDB ID: 2XQ2) [41] , isolated from Vibrio parahaemolyticus, has the highest identities and similarities with hSGLT2 sequences, it was selected as the template to build the homology model of hSGLT2 using Build Homology Models of DS 2017 R2. Subsequently, the homology model obtained after 10 ns of MD refinement was evaluated via Ramachandran plot [42] on RAMPAGE website and Verify Protein (Profiles-3D) protocol of DS 2017 R2 to assess the compatibility of the 3D structure with the corresponding amino acid sequence [43] .

2.6. Molecular Docking

Molecular docking was used to evaluate the binding affinity and clarify the binding mode of the screened compounds. Since vSGLT (PDB ID: 3DH4) [44] is the only co-crystallized structure of vSGLT and galactose, the position of galactose was used to identify the binding site of hSGLT2, with a sphere radius of 10 Å around the binding site being defined as the docking region [45] . CDOCKER program [46] of DS 2017 R2 with the CHARMm force field was used for precise and flexible docking. All of the docking parameters were set to default. For each screened compound, 10 docking poses were generated and the pose with the lowest CDOCKER interaction energy was compared to that of canagliflozin. The compounds with lower CDOCKER interaction energy than that of canagliflozin were finally selected.

2.7. ADMET Prediction

In this study, the online sever admetSAR was employed to predict the ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties of the screened compounds. It contains 22 qualitative classification and 5 quantitative regression models with highly predictive accuracy, allowing to estimate ADMET properties of the inputted compounds [47] . Each selected compound was converted to the simplified molecular-input line-entry systems (SMILES) format and inputted into the admetSAR for prediction. The compounds with superior ADMET properties were considered as potential hSGLT2 inhibitors.

2.8. Molecular Dynamics (MD) Simulation

The binding stabilities between hSGLT2 structure and the selected compounds were evaluated by MD simulations using Gromacs 2016.3 [48] . Combined with NVIDIA CUDA 8.0, Antechamber in AMBER tool 2017 and ACPYPE, the AMBER99SB-ILDN force field was applied for the hSGLT2 and the general AMBER force field (GAFF) was applied for the screened compounds in the MD simulations process [49] [50] . All systems were immersed in a 10.0 Å dodecahedron periodic box with TIP3P water molecules. The sodium and chloride ions were added to reach 0.15 M as physiology state. The temperature was heated gradually from 0 to 310 K and then equilibrated to 310 K at 1.0 atm pressure. The entire system was stabilized through energy minimization and then executed two steps of restrained equilibrations for 200 ps, including NVT equilibration and NPT equilibration. Subsequently, the full MD simulations were performed for 150 ns in the equilibrate state.

2.9. Binding Free Energy Calculation

The Poisson-Boltzmann and surface area continuum solvation (MM/PBSA) method implemented in g_mmpbsa for GROMACS-5.1.X was used to calculate the binding free energies (ΔG_binding) between hSGLT2 and the selected compounds [51] . ΔG_binding was composed of four individual components, including van der Waals (ΔG_vdw), electrostatic (ΔG_elec), polar (ΔG_polar) and nonpolar (ΔG_nonpolar) solvation energy. The solvent accessible volume only (SAV-only) nonpolar model was applied for nonpolar solvation energy calculation in MM/PBSA. The temperature was set to 310 K, while the other parameters were set to default. For a more precise binding free energy analysis, the stable last 20 ns in the MD simulation trajectory was selected to generate a total of 2000 snapshots.

3. Results and Discussion

In the early stage of hSGLT2 inhibitor development, the strategies were committed to chemical synthesis, structure-activity relationships (SAR) analysis, and lead modification. However, no significant improvements in efficacy have been made due to the highly structure similarity. Through the aid of several computational approaches, precise molecular level study of hSGKT2-inhibitors binding interaction has become possible. To the best of our knowledge, this study is the first attempt to discover novel N-glycoside and non-glycoside hSGLT2 inhibitors through a combination of several computational approaches, including ligand-based pharmacophore modeling, virtual screening, homology modeling, molecular docking, ADMET predictions, molecular dynamics (MD) simulations and binding free energy calculations.

3.1. Pharmacophore Model Generation

Pharmacophore modeling has been proven to be very successful not only in demonstrating the structure-activity relationships between ligands and receptors, but also in developing new drugs. A pharmacophore model can be established either in a ligand-based manner or in a structure-based manner. While structure-based pharmacophore model works directly with the 3D structure of the target or its complex with a ligand, ligand-based pharmacophore model is derived from the structure of a series of ligands by analyzing the common pivotal features of the superimposed ligands and translated into pharmacophore features subsequently. The generated pharmacophore model offers a rational hypothetical view of these significant chemical features responsible for activity. Since the crystal structure of hSGLT2 is not yet available, direct structure-based pharmacophore modeling is not applicable.

In this study, a total of 10 pharmacophore hypotheses were established based on the training set compounds (Figure 1) using the HypoGen module embedded in DS 2017 R2 and the results are given in Table 1. To choose the best hypothesis, a cost analysis was assessed on the basis of two important theoretical cost calculations (represented in bit units): the fixed cost, which represents the simplest model that fits all data perfectly, and the null cost, which represents the highest cost of a model with no correlation. The total cost of an excellent pharmacophore model should be close to the fixed cost and far away from the null cost. The RMSD value shows the quality of the correlation between the experimental and the estimated activity data. To verify the predictive ability of all hypotheses on the training set compounds, the activity of each compound was estimated by the regression analysis [52] . According to the above mentioned criteria, the first hypothesis (Hypo1) was chosen to be the best hypothesis characterized by the lowest total cost (110.381), the highest correlation coefficient (R = 0.919) and the lowest RMSD value (0.993).

Figure 3(A) and Figure 3(B) show the 3D spatial arrangement and distance constraints of the features (three HBDs and one HY) in Hypo1. The features of Hypo1 were mapped onto the most active compound 1 with a fit value of 10.616 and the least active compound 22 with a fit value of 7.452 as shown in Figure 3(C) and Figure 3(D), respectively. Fit value indicates how well the features in a pharmacophore model overlap the chemical features in the molecule and thus

^aCost difference = null cost − total cost. null cost = 140.937. Fixed cost = 97.198. All cost values are in bits; ^bHBA, hydrogen bond acceptor; HBD, hydrogen bond donor; HY, hydrophobic; RA, ring aromatic; EV, excluded volumes.

aid in understanding the chemical meaning of the hypothesis. The difference between the fit values of these two compounds may explain the difference in their activities.

3.2. Pharmacophore Model Validations

The accuracy of Hypo1 is important to achieve success in the subsequent virtual screening process, therefore its validity was assessed by test set prediction and GH score method.

3.2.1. Test Set Prediction

The significance of the selected pharmacophore model depends on its ability to predict the biological activity of test set compounds along with the training set compounds. A total of 23 structurally diverse test set compounds with a wide range of activities were utilized to verify the prediction ability of Hypo1. All the test compounds were subjected to the same conformation generation procedures as for the training set compounds. Table 2 shows the experimental and estimated activities of the 22 training set compounds and Table 3 shows the experimental and estimated activities of the 23 test set compounds.

To verify the predictability power of Hypo1, the estimated activity (pIC₅₀) for each individual test set compound was calculated in order to correlate to its experimental activity using simple regression analysis. The obtained correlation

^aError factor calculated as the ratio of the experimental activity to the estimated activity; “+” indicates that the estimated IC₅₀ is higher than the experimental IC₅₀; “-” indicates that the estimated IC₅₀ is lower than the experimental IC₅₀; ^bFit value indicates how well the compound mapped onto the features of Hypo 1; ^cActivity scale: highly active (IC₅₀ < 10 nM, +++), moderately active (10 nM ≤ IC₅₀ < 1000 nM, ++), and inactive (IC₅₀ ≥ 1000 nM, +).

^aError factor calculated as the ratio of the experimental activity to the estimated activity; “+” indicates that the estimated IC₅₀ is higher than the experimental IC₅₀; “-” indicates that the estimated IC₅₀ is lower than the experimental IC₅₀; ^bFit value indicates how well the compound mapped onto the features of Hypo 1; ^cActivity scale: highly active (IC₅₀ < 10 nM, +++), moderately active (10 nM ≤ IC₅₀ < 1000 nM, ++), and inactive (IC₅₀ ≥ 1000 nM, +).

coefficient values were 0.962 and 0.919 for the test set and training set compounds, respectively (Figure 4). The results demonstrate that Hypo1 is a reliable pharmacophore model with high predictability power.

3.2.2. GH Score Method

The main purpose of the GH score method is to assess the discriminative ability of Hypo1 to correctly identify the known active compounds from the inactive ones. All the detailed GH score parameters of Hypo1 are summarized in Table 4. To perform GH score calculation, an external database, consisting of a total of 550 compounds (D) including 50 actives (A), was utilized to evaluate the discriminative ability of Hypo1. After screening this database by Hypo1, 57 compounds were retrieved as hits. Among them, 49 compounds were from the known actives (Ha). According to the previous study, the GH score of a good screening model must be higher than 0.7 [53] . The GH score of Hypo1 was calculated to be 0.876, indicating that it has the ability in selecting the active compounds among inactive ones.

3.3. Virtual Screening

A schematic flowchart of the overall virtual screening process in this study is provided in Figure 2. To effectively screen the NCI and TCM databases with 265,242 and 60,563 compounds, respectively, Lipinski rule of five [39] and Veber

^aD is the number of total compounds in database, including 50 active compounds and 500 decoys; ^bHt is the total number of hit compounds in database, including 49 active compounds and 8 decoys; ^{c $\left[\left(\frac{\text{Ha}}{4\text{Ht}\times \text{A}}\right)\left(3\text{A}+\text{Ht}\right)\times \left(1-\left(\frac{\text{Ht}-\text{Ha}}{\text{D}-\text{A}}\right)\right)\right]$} , GH score > 0.7 indicates a good pharmacophore model.

rule [40] were first employed to ensure that the screened compounds are drug-like. The remaining 214,792 and 19,026 compounds that passed these rules were subsequently mapped onto Hypo1 to identify potential hSGLT2 inhibitors. Finally, 1617 and 1193 compounds from NCI and TCM databases, respectively, which pass all the criteria mentioned above were subjected to further molecular docking studies. In this stage, our Hypo1 model shows great capability of efficiently identifying potential hits from a large number of compounds.

3.4. Homology Modeling

Homology modeling is an important method of generating the 3D structure of a target protein from a closely related template. The accuracy of the obtained structural model dramatically depends on the identity and similarity between the target and the template proteins. Our sequence alignment results indicate that the template vSGLT (PDB ID: 2XQ2) has a sequence identity of 24.0% and a similarity of 47.0% with hSGLT2. It is generally assumed that the sum of the BLAST sequence identity and similarity of at least 60% is suggestive of an effective sequence template. Thus, the crystal structure of vSGLT was herein adopted as the template and the pair-wise sequence alignment was used to construct the homology model of hSGLT2. The 3D structural model of hSGLT2, including 643 amino acids, was generated using Build Homology Models of DS 2017 R2. The model obtained after 10 ns of MD simulation refinement was further validated using a Ramachandran plot (RC) and the Verify Protein (Profiles-3D) protocol. The corresponding RC plot showed that 99.5% of the residues in the hSGLT2 structural model were found in both the most favored and additional allowed regions, and only 0.5% of the disfavored residues were found to be located in the loop part of the protein. It indicates that the torsion angles of the residues in the homology model are reliable [54] . The verification score obtained from Profile-3D upon inspecting the quality of the protein structure geometry was 193.37, close to the expected high verification score of 233.98, indicating that this homology model yields satisfactory results. We further examined that some structural features of this homology model were consistent with the characteristics of the SGLT family, which includes 14 transmembrane helices with inverted repeat topology segments (Figure 5).

3.5. Molecular Docking

Molecular docking has been intensively used to clarify the binding mode of the compounds and to obtain the other information that could be utilized for further structural optimization. Although the structural model of hSGLT2 was built by homology modeling successfully, the binding site of hSGLT2 is still undefined. So far, vSGLT (PDB ID: 3DH4) [44] is the only co-crystallized structure of vSGLT and galactose. Therefore, the position of galactose was used to identify the binding site of hSGLT2 in this study. To distinguish the binding affinity of each screened compound, CDOCKER interaction energy (kcal/mol) was used as the selecting criteria. Three approved hSGLT2 inhibitors, including canagliflozin, dapagliflozin and empagliflozin, were used to determine the possible binding modes of the selected hSGLT2 inhibitor. These three hSGLT2 inhibitors show the identical binding mode in the binding site of hSGLT2 (Figure 6(A)). Among them, canagliflozin has the lowest CDOCKER interaction energy of −47.66 kcal/mol (Figure 6(C)). Then, the 1607 and 1193 compounds selected from NCI and TCM databases, respectively, were docked into the binding site of hSGLT2 and the compounds with lower CDOCKER interaction energies than

that of canagliflozin were selected as potential hits, resulting in 332 and 175 compounds from NCI and TCM databases, respectively. These selected compounds were subjected to visual inspection to confirm their binding modes and whether their structural features belong to either N-glycoside or non-glycoside.

Previous studies have indicated that the hSGLT2 inhibitors fit into the binding site of hSGLT2 by interacting with several important residues, such as Asn75, His80, Glu99, Ala102, Arg267, Tyr290, Trp291, Gln457, and Pro518 [55] [56] [57] [58] . Among these important residues, Asn75, Glu99, Arg267, and Tyr290 are the most essential key residues since they form a hydrogen bonding network with hSGLT2 inhibitors [55] [56] [57] [58] . According to the above screening criteria, only 1 N-glycoside and 4 non-glycoside compounds were finally identified to be potential hSGLT2 inhibitors (Table 5). These screened compounds combined with essential key residues and formed stable combinations through hydrogen bonds. These selected compounds, which exhibit similar binding mode and higher binding affinities than canagliflozin, were further evaluated through ADMET prediction followed by MD simulation studies.

3.6. ADMET Prediction

Prediction of the compound ADMET properties describes the disposition of a pharmaceutical compound within the human body. Considering that most of the previously reported hSGLT2 inhibitors failed in their clinical trials due to poor ADMET properties, filtering and optimization of the ADMET properties in the early stage of the drug development are necessary to avoid failure. In this study, the ADMET properties of canagliflozin and the selected 1 N-glycoside (NSC679207) and 4 non-glycoside (NSC727709, NSC231775, TCM21610 and TCM_Piperenol_A) compounds were predicted by the online sever admetSAR and the results are presented in Table 6. In terms of absorption, only NSC231775 have bad property in human intestinal absorption (HIA), suggesting oral delivery is infeasible. In terms of toxicity, canagliflozin, NSC679207, TCM21610 and TCM_Piperenol_A are all non-mutagenic (AMES toxicity test) and non-carcinogenic (carcinogenic profile) [59] . In human heart, hERG potassium channels are essential for regulating electrical activity. The blockage of hERG channel by drug causes lethal long QT syndrome (LQTS), which leads to withdraw from the clinical trials or markets [60] . The results of hERG inhibitor prediction I/II showed that TCM21610 is an hERG inhibitor, involving high LQTS risk than the others.

Since the US Food and Drug Administration (FDA) published the first in vitro/in vivo drug interaction guidance documents in 1997 and 1999, respectively, off-target transporter interactions and drug-drug interactions (DDIs) have become the major challenges for drug development. CYP450 enzymes-based unanticipated DDIs and drug metabolism problems are also a common cause of adverse drug events (ADE). Interacting with CYP450 enzymes may lead to DDIs, hepatotoxicity, reactive metabolite formation, idiosyncratic adverse drug interactions, and/or even loss the efficacy of drugs. CYP1A2 is related to alcohol metabolism, indicating that the inhibition of CYP1A2 results in the reducing of alcohol clearance and increases the risk of liver toxicity [61] . The inhibition of CYP2C19 will affect the metabolism of other drugs through the CYP2C19 pathway,

^a+ for positive and - for negative; ^b+ for substrate and - for non-substrate; ^c+ for inhibitor and - for non-inhibitor; ^d+ for penetrable and - for non-penetrable; ^e+ for toxic and - for non-toxic; ^f+ for carcinogens and - for non-carcinogens.

which is related to the function of platelet in human body. Therefore, CYP2C19 inhibitors should be carefully used while cooperating with drugs that affect platelet aggregation [62] . The level of CYP inhibitory promiscuity inhibition represents the risks of causing CYP450 enzymes-based unanticipated DDIs and metabolism problems [63] . According to the above information, canagliflozin and NSC727709 must pay attention to unanticipated DDIs and the risk of hepatotoxicity.

Based on the ADMET prediction, among the five selected compounds, only NSC679207 and TCM_Piperenol_A exhibit better ADMET properties than that of canagliflozin, indicating that both compounds exhibit great potential to serve as N-glycoside and non-glycoside hSGLT2 inhibitors.

3.7. MD Simulations

MD simulations are able to analyze the trajectory of protein-ligand complex during the dynamics process to validate the stability of docking results and calculate binding free energies in the entire process. The structures of hSGLT2 docked with canagliflozin, NSC679207, and TCM_Piperenol_A were further subjected to 150 ns MD simulations to investigate the binding stabilities of these three compounds towards hSGLT2. The backbone RMSD of the three compounds and hSGLT2 revealed that an equilibrated and converged state was achieved after 50 ns of simulation. Therefore, MD simulations were then extended to 150 ns in an effort to gauge the stabilities of these three binding systems [64] . As shown in Figure 7(A), no major fluctuations were observed. Canagliflozin, NSC679207, and TCM_Piperenol_A were maintained at stable state with averaged RMSD values of 3.195 Å, 3.178 Å, and 2.973 Å, respectively in the last 10 ns MD simulations. It allows the changes in the conformations of these three compounds and the residues of the binding site to be discerned.

To explore the protein residues flexibility throughout the MD simulation, the root-mean-squared fluctuations (RMSF) of individual hSGLT2 residues were calculated. As shown in Figure 7(B), NSC679207 and TCM_Piperenol_A exhibited similar RMSF distributions comparing to canagliflozin. In addition, the conformational changes of the above mentioned hSGLT2 key residues (Asn75, Glu99, Arg267, and Tyr290) were all very small (Figure 7(C)). Our results show that both NSC679207 and TCM_Piperenol_A exhibit stable RMSF distribution, indicating that they can sufficiently fix most of the hSGLT2 key residues (Figure 7(C)). The results of RMSF analyses suggest that both NSC679207 and TCM_Piperenol_A have extremely stable binding stabilities, similar to that of canagliflozin, indicating

that they form tight binding complexes with hSGLT2 during the entire MD simulation processes.

3.8. Binding Free Energy Calculation and the Analysis of Binding Interactions

In the past, molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) or molecular mechanics/generalized Born surface area (MM/GBSA) methods have been intensively applied to determine the correlations between the calculated free energy and the experimental activity. In this study, the binding free energy was calculated using the MM/PBSA method implemented in GROMACS-5.1.X 2000 snapshots from the last 20 ns of the MD simulations of the complexes of hSGLT2 with canagliflozin, NSC679207, and TCM_Piperenol_A were used to perform binding free energy (ΔG_binding) calculations and the results are as summarized in Table 7. NSC679207 (−45.30 kcal/mol) and TCM_Piperenol_A (−55.19 kcal/mol) showed higher and lower binding free energies, respectively, comparing to that of canagliflozin (−52.32 kcal/mol).

All glycoside hSGLT2 inhibitors contain a glycoside structure, which binds to the glucose binding site of hSGLT2 and an aglycone tail, which binds to the wall of the hydrophilic cavity leading to enhance the stability [65] . The most important binding affinities of the hSGLT2 inhibitors were provided by hydrogen bonding network around the glycoside structure. To understand the difference of the binding free energies among canagliflozin, NSC679207, and TCM_Piperenol_A, the initial binding interactions of these three compounds were revealed through their 3D structures and 2D binding interactions diagrams (Figure 8). According to Figure 8(A), Figure 8(C), and Figure 8(E), these three compounds exhibit the identical binding mode and stably locate in the binding pocket of hSGLT2 through the above mentioned hydrogen bonding network. The 2D interactions (Figure 8(B), Figure 8(D), and Figure 8(F)) also showed that these three compounds are all able to form the essential hydrogen bonds with hSGLT2 through several key residues, including Asn75, Glu99, Arg267, and Tyr290. The difference of the binding free energies between canagliflozin and NSC679207 (N-glycoside compound) are mostly occurred in the term of ΔG_vdw, indicating that the difference may be mainly attributed in the aglycone side chain instead of the competitive binding site. It indicates that although NSC679207 has higher binding free energy than that of canagliflozin, it is still sufficient to serve as a novel N-glycoside hSGLT2 inhibitor. On the other hand, TCM_Piperenol_A (non-glycoside compound) showed great potential to be a promising non-glycoside hSGLT2 inhibitor

due to its good ADMET preperties, binding stability, and better binding free energy than canagliflozin.

To the best of our knowledge, this is the first attempt to employ a combination of a series of computational approaches to discovery novel N-glycoside and non-glycoside hSGLT2 inhibitors. In this study, a ligand-based pharmacophore containing 3 HBDs and 1 HY features was constructed and validated. The subsequent docking results showed that the most impotrtant binding affinities of the selected hSGLT2 inhibitors were provided by hydrogen bonding network, which is consistent with the previous studies. TCM is by far the largest free downloadable database for Chinese herbal medicine in the world. In previous studies, the possibility of natural compounds to serve as hSGLT2 inhibitors have been revealed ever since the first SGLTs inhibitor, phlorizin, a substance found in some fruit trees, was discovered [66] . TCM_Piperenol_A is isolated from Piper cubeba (cubeb pepper), which is mainly grown in the islands of Java and Sumatra and usually used as a spicy ingradient. Our study establishes a fast and low cost platform to successfully discover novel and N-glycoside and non-glycoside hSGLT2 inhibitors from various datebases, which may become potential leads for the treatment of T2DM after a series of in vitro and/or in vivo bioassays in the future.

4. Conclusion

The priority of discovering novel hSGLT2 inhibitors with less side effects has been emphasized and increased year by year for the development of T2DM drugs [67] . However, the development of hSGLT2 inhibitors is still limited to the C-glycoside inhibitors in both commercially available drugs and clinical trials. Although the advantages of N-glycoside and non-glycoside inhibitors have been reported, the discovery of novel drugs is still high cost and time-consuming. In this study, we used highly accurate pharmacophore model and virtual screening to retrieve various types of compounds from large databases rapidly. In addition to the well established NCI database, the TCM database, which contains abundant Chinese herbal medicines and natural products, was also used as the screened database in this study. Finally, NSC679207 and TCM_Piperenol_A, which exhibit different chemical scaffolds than that of canagliflozin, were successfully identified to possess higher binding affinities towards hSGLT2. Besides, they also exhibit good pharmacokinetic properties and low risk in several common lethal side effects via ADMET prediction. Thus, these two compounds exhibit great potential to serve as novel N-glycoside and non-glycoside hSGLT2 inhibitors and may be further subjected to a series of in vitro and/or in vivo bioassays. Finally, our results also open a new channel leading to a prospective development for T2DM therapy.

Acknowledgements

The authors greatly thank the Ministry of Science and Technology of Taiwan (MOST 106-2221-E-027-117-MY3 and MOST 107-2221-E-027-038-MY3) and National Taipei University of Technology and Taipei Medical University joint project (USTP-NTUT-TMU-105-05) for their financial supports.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Chang, C.-Y., Ho, Y., Lin, S.-J. and Liu, H.-L. (2019) Discovery of Novel N-Glycoside and Non-Glycoside hSGLT2 Inhibitors for the Treatment of Type 2 Diabetes Mellitus. Journal of Diabetes Mellitus, 9, 77-104. https://doi.org/10.4236/jdm.2019.93009

References