In-silico screening and identification of phytochemicals from Centella asiatica as potential inhibitors of sodium-glucose co-transporter 2 for treating diabetes

Abstract

Sodium-glucose co-transporter 2 (SGLT-2) is a major transport protein responsible for reabsorption of glucose from the kidney back to the bloodstream. Inhibiting this protein effectively lowers the glucose level of diabetic patients; however, the use of synthetic SGLT-2 inhibitors has been linked to some serious adverse effects. There is a need to identify safer alternatives that are equally or more effective as the current inhibitor drugs. Phytochemicals are known for their efficacy as herbal remedies, but these molecules remain underexplored as source of therapeutic agents. In this study, we performed in silico screening to identify potential SGLT-2 inhibitors from the 21 phytochemicals from Centella asiatica. Docking results identified eleven compounds with estimated binding energies comparable to that of known inhibitors drugs. The stability of the complexes was then elucidated using 100 ns MD simulations. From our dynamic binding free energy calculations using MM/PBSA, asiaticoside, betulinic acid, centellasapogenol, methyl brahmate, and rutin exceeded at least one of the binding energies of the reference compounds, which highlights their strong affinity towards SGLT-2. Among the five, betulinic acid, centellasapogenol, and methyl brahmate maintained their structural stability to the same extent as the references and exhibited better oral bioavailability and excellent drug-like properties. Because of these results, it is recommended to prioritize betulinic acid, centellasapogenol, and methyl brahmate in future in vitro and in vivo studies to verify their potential as inhibitor drugs for diabetes therapies.

Communicated by Ramaswamy H. Sarma

Keywords:

• Centella asiatica
• diabetes
• SGLT-2
• molecular dynamics
• phytochemicals
• molecular docking
• ADMET

Introduction

Diabetes is a chronic disease characterized by an increased level of glucose in the blood due to the body’s inefficiency to produce or use insulin. It is mostly common to adults with 8.5% prevalence and tends to be diagnosed in older people in their 40–50s (Centers for Disease Control and Prevention, 2020). If left uncontrolled, diabetes can lead to several complications including heart attack, kidney disease, lower-limb amputations, and adult-onset blindness. Additionally, 43% of diabetic patients die prematurely before the age of 70 due to complications (Centers for Disease Control and Prevention, 2020; Roglic, 2016). In the Philippines, diabetes has one of the highest mortality rates, accounting for 33,295 deaths in 2016, and it is also the 6^th leading cause of death in the country (Bersales, 2018).

Sodium-glucose co-transporter 2 (SGLT-2) inhibitors are among the newest class of anti-hyperglycemic medications to receive approval from the FDA (Hsia et al., 2017). SGLT-2 is the most prevalent sodium-glucose co-transporter subtype which facilitates 90% of glucose reabsorption in the kidney back to the bloodstream (Hsia et al., 2017). It is therefore an ideal target for the treatment of diabetes because of their large contribution to glucose reabsorption. SGLT-2 inhibitors have a unique mechanism that is independent of the level of insulin secretion–they block glucose attachment to the protein in the kidney and disrupt the main function of the sugar transporter thereby preventing glucose reabsorption to the bloodstream. Once SGLT-2 transport of glucose is inhibited, glucosuria or the excretion of glucose through the urine is induced to lower the sugar level of the body.

The binding site of the SGLT-2 protein is sandwiched by a group of residues that form intracellular and extracellular gates. In the intracellular gate, Tyr 263 forms the strongest van der Waals interactions to hold the glucose molecule and prevents the sugar’s exit further into the cavity (Tamura et al., 2015). The extracellular gate, composed of Met 73, Tyr 87, and Phe 424, hinders the detachment of glucose (Wright et al., 2011). Certain residues such as Ala 294, Gln 428, Gln 69, Glu 68, Glu 88, Ser 91, and Asn 260 were also determined to be of importance in the absorption mechanism (Tamura et al., 2015; Wright et al., 2011). Among these, Gln 428 is established to be the most vital because of its direct interaction with the sugar for absorption and any mutation in this residue hinders sugar absorption (Rizvi et al., 2013). Additionally, Glu 88 and Glu 68 provides the strongest electrostatic interactions – a significant contributor to the strong binding affinity observed between the sugar molecule and the protein (Tamura et al., 2015)

There are currently three approved SGLT-2 drugs by the FDA – canagliflozin (CAN), dapagliflozin (DAP), and empagliflozin (EMP) – while some are still undergoing clinal trials (Hsia et al., 2017). However, the use of these drugs are associated to the development of urinary tract infections and genital infection (Singh & Kumar, 2018). Aside from these common adverse effects of SGLT-2 inhibitors, the chemically synthesized SGLT-2 inhibitor drugs carry additional health risks. An increased risk of bone fractures and decreased bone mineral density have been linked to the use of canagliflozin, for which the FDA issued a relevant warning (Tentolouris et al., 2019). Lower limb amputations and pancreatitis – though exceedingly rare – were also observed from the use of this drug. For dapagliflozin, patients reportedly developed nasopharyngitis, a unique adverse reaction of the drug (Sha et al., 2015), while empagliflozin has been correlated with an increased chance of upper respiratory infections which led to a drug label warning from the FDA. The exact mechanism of this adverse effects remains unknown so whether this is caused by SLGT-2 inhibition or not is yet to be verified (Tentolouris et al., 2019). Because of the occurrence of these undesirable side effects, there is a need for novel agents that may offer superior efficacy and safer treatment options.

One alternative source of drug candidates are phytochemicals or bioactive compounds from plants whose vast structural diversity of natural molecules remains largely untapped. 80% of the population from developing countries still relies on herbal extracts because of their proven safety and efficacy in traditional medicine (Veeresham, 2012). Natural compounds have also been established in various in vivo and in vitro studies as one of the major sources of anti-diabetic agents (Jugran et al., 2020; Kabir et al., 2014; Supkamonseni et al., 2014). The first SGLT-2 inhibitor discovered, phlorizin, was a natural compound isolated from the bark of an apple tree; however, the work was discontinued because of the compound’s non-selective nature and poor oral bioavailability (Kalra, 2014). Following this, numerous studies began characterizing extracts from plants as anti-diabetic inhibitors or activator of various proteins but only few explored the inhibition of the SGLT-2 protein (Choi, 2016; Hung et al., 2012). Thus, efforts to identify novel SGLT-2 inhibitors derived from plants are still on-going.

In this study, we explored the potential of the phytochemicals from Centella asiatica as an SGLT-2 inhibitor (Figure 1). Centella asiatica is a small and edible perennial plant with a wide range of folk medicinal applications indigenous to Asian and African countries. Extracts from the plant have been found to exhibit antioxidant (Pittella et al., 2009), antifungal, antibacterial (Dash et al., 2011), and anti-inflammatory (George & Joseph, 2009) properties mostly attributed to its vast number of phytochemicals, specifically triterpenes and flavonoids (Hashim et al., 2011). Although Centella asiatica’s main usage in folk medicine is to treat skin lesions and promote wound healing (Gohil et al., 2010), some in-vivo studies have probed the use of the plant’s extract as a cure for diabetes. A significant reduction of glucose by up to 69% was observed in induced diabetic rats from the methanolic extract of the plant which was highly comparable to the effects of the known diabetic drug, glibenclamide (Chauhan et al., 2010). There was also a protective effect against diabetes which may be due to the plant’s high number of flavonoids working synergistically. Two more recent studies identified the triterpene, asiaticoside (AST), and the flavonoid, rutin (RUT) as the active ingredients responsible for the anti-diabetic activity observed in rats. It was also theorized that asiaticoside enhanced the insulin secretion of pancreatic β-cells that led to the reduced glucose levels (Fitrianda et al., 2017). Rutin, on the other hand, was shown to slow down carbohydrate digestion and glucose absorption in the body (Supkamonseni et al., 2014). But the exact pharmacological mechanisms of these effects are not yet clear, thus glycemic control due to SGLT-2 inhibition could be a possibility. The continuous use of Centella asiatica in herbal medicine, and its perceived therapeutic efficacy and versatility has driven the scientific community to accumulate large bodies of evidence to highlight the immense potential of this plant to cure various diseases (Gohil et al., 2010). To the best of our knowledge, an investigation of the inhibitory properties of Centella asiatica’s phytochemicals against SGLT-2 have never been done. This study aims to test if Centella asiatica is a potential source of safer alternative drug molecules capable of inhibiting SGLT-2 to lower the chances of diabetes complications.

Figure 1. Two-dimensional structures of the 21 phytochemicals of Centella asiatica considered in this study. Also shown are the structures of canagliflozin, dapagliflozin, and empagliflozin – the known SGLT-2 inhibitors.

Materials and methods

The phytochemical structures of triterpenes and flavonoids (Table S1, Supplementary Materials) of Centella asiatica were obtained from PubChem (www.pubchem.com) (Kim et al., 2016). The structures of brahmol (BRA) and methyl brahmate (MEB) are unavailable in PubChem so their structures were constructed in Maestro (Schrödinger, 2020a) and were prepared in Avogadro 1.2 (Hanwell et al., 2012) using General AMBER Force Field (GAFF). Steepest descent algorithm with 10,000 iterations was employed to optimize and establish the relaxed conformation. In total, 21 triterpenes and flavonoids were considered in this study. The protein structure of SGLT-2 (PDB ID: 2XQ2) (Watanabe et al., 2010) (Figure S1, Supplementary Materials) was imported from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (www.rcsb.org) (Rose et al., 2016).

Table 1. ADMET and drug-likeness results of the ligands and the references.

Ligands	Lipinski’s rule of five					Selected QikProp descriptors							Molsoft drug-likeness score	Toxicity
	Mol MW	Donorhb	accpthb	QPlogPo/w	No. of violations	HOA	PHOA	QPlogS	QPlogPC16	QPlogPw	QPlogPoct	QikProp score
Reference range	<500	≤5	≤10	<5		1: low, 2: medium, 3: high	>80%: high <25%: poor	–6.5 − 0.5	4.0 − 18.0	4.0 − 45.0	8.0 − 35.0	max: 45	0–2
Canagliflozin	444.517	4	8.5	3.144	0	3	88.319	–5.396	14.652	17.316	26.11	41	0.28	PASS
Dapagliflozin	408.878	4	9.25	2.088	0	2	81.349	–4.389	13.674	16.948	24.08	41	0.48	PASS
Empagliflozin	450.915	4	10.95	1.922	0	2	80.374	–4.391	14.34	18.625	26.398	41	0.62	PASS
Apigenin	270.241	2	3.75	3.75	0	3	73.64	–3.374	9.716	2	14.459	41	0.77	PASS
Asiatic acid	488.706	4	7.1	4.173	0	3	84.759	–5.148	13.361	14.311	25.622	42	0.71	PASS
Asiaticoside	959.133	12	30.9	–1.539	3	1	0	–3.547	27.417	48.87	63.15	23	0.64	PASS
Bayogenin	488.706	4	7.1	4.27	0	3	85.457	–5.363	13.498	14.343	25.902	43	0.59	PASS
Betulinic acid	456.707	2	3.7	6.2	1	1	95.645	–6.712	12.29	8.061	20.741	39	0.31	PASS
Brahmol	490.722	5	8.5	3.103	0	3	92.263	–4.776	13.553	16.7	27.265	41	0.07	PASS
Castilliferol	432.386	3	7	2.618	0	2	68.689	–5.354	15.562	15.672	23.323	40	0.95	PASS
Centellasapogenol	492.738	5	8.5	3.099	0	3	89.839	–4.462	13.603	16.606	27.359	42	0.47	PASS
Isothankunic acid	504.706	5	7.85	3.879	1	3	69.864	–4.887	13.803	16.302	27.776	42	0.91	PASS
Kaempferol	286.24	3	4.5	1.053	0	3	64.189	–3.12	10.216	12.315	16.613	41	0.77	PASS
Madasiatic acid	488.706	4	7.1	4.338	0	3	82.124	–5.829	13.901	14.678	26.942	43	0.53	PASS
Madecassic acid	504.706	5	8.8	3.255	1	2	60.838	–4.672	13.952	17.481	28.615	40	0.69	PASS
Madecassoside	975.132	13	32.6	–2.798	3	1	0	–2.108	26.767	51.985	64.923	25	0.62	PASS
Methyl asiatate	502.733	3	7.1	4.684	1	3	94.263	–6.132	13.458	12.586	25.021	44	0.58	PASS
Methyl brahmate	532.759	3	8.8	4.438	1	3	95.836	–6.142	13.897	14.093	26.131	43	0.72	PASS
Myricetin	318.239	5	6	–0.276	1	2	27.77	–2.626	11.219	16.484	20.664	40	–0.04	PASS
Patuletin	332.266	4	6	0.628	0	2	57.292	–3.007	10.97	14.593	19.068	40	0.78	PASS
Pomolic acid	472.707	3	4.45	5.637	1	3	90.982	–6.13	12.513	10.228	22.578	41	0.29	PASS
Quercetin	302.24	4	5.25	0.362	0	2	52.349	–2.83	10.666	14.379	18.492	40	0.93	PASS
Rutin	610.524	9	20.55	–2.638	3	1	0	–1.785	18.235	35.602	41.721	30	1.1	PASS
Terminolic acid	332.266	4	6	0.628	0	2	57.292	–3.007	10.97	14.593	19.068	40	0.54	PASS

The Qikprop descriptors from left to right: molecular weight (mol MW), H-bond donors (donorHB), H-bond acceptors (accptHB), predicted octanol/water partition coefficient (QPlogPo/w), human oral absorption (HOA), percentage human oral absorption (PHOA), predicted aqueous solubility (QPlogS), predicted hexadecane/gas partition coefficient (QPlogPC16), predicted octanol/gas predicted water/gas partition coefficient (QPlogPw), and predicted octanol/gas partition coefficient (QPlogPoct).

ADMET and drug-like property evaluation

The QikProp module of Schrödinger (Schrödinger, 2020b) was used to evaluate the absorption, distribution, metabolism, elimination, and toxicity (ADMET) properties of all the 21 phytochemicals and the three known SGLT-2 inhibitors. There are 45 descriptors and parameters from the QikProp modules, each with their own range of accepted values. Part of the QikProp descriptors is Lipinski’s Rule of Five (Ro5), a simple yet widely recognized rule of thumb to evaluate the drug’s oral absorption. The rule states that an ideal candidate drug must have less than five H-bond donors and less than ten H-bond acceptors, must weigh less than 500 g/mol, and the calculated octanol-water partition coefficient log P must be less than five. Having more violations from these guidelines increases the likelihood of poor solubility and permeability (Lipinski, 2004). Two additional servers were also used to evaluate the compounds’ drug-likeness and toxicity. Molsoft (www.molsoft.com) (Molsoft, 2018) provides the overall drug-likeness score and ChemBioSever (www.chemibioserver.com) (Athanasiadis et al., 2012) detects the presence of some known toxic moieties.

Molecular docking

Molecular docking was performed in Autodock 4.2 (Morris et al., 2009) using Lamarckian Genetic Algorithm (LGA). The ligands were treated as flexible molecules and the receptor was kept rigid. The grid center coordinates for the binding cavity of SGLT-2 were set to 2.111, 37.221, 52.376 in the x, y, z direction, and a grid box size of 40 × 40 × 40 Å with 0.375 Å spacing was used to cover the known active residues within the binding cavity. The location of the active site of SGLT-2 was determined using the known residues that surround the active pocket of the protein (Deng et al., 2014; Wright et al., 2011). For exhaustiveness, the number of GA runs was set to 1000, with 1000 poses generated per run. Only the best pose per run was saved. This is to get as close as possible to the correct conformation. The docking score of each pose was estimated using the binding free energy (ΔG, Equation (1)(1) $$\mathrm{ }\Delta G=\Delta G_{\mathit{vdW}}+\Delta G_{\mathit{hbond}}+\Delta G_{\mathit{desolv}}+\Delta G_{\mathit{elec}}+\Delta G_{\mathit{tor}}$$(1) ), which is the sum of the energies from van der Waals and H-bond interactions, desolvation energy, electrostatic energy, and torsional free energy (Morris et al., 2009). Aside from the 21 phytochemicals, three known SGLT-2 inhibitors are used as control compounds. (1) $$\mathrm{ }\Delta G=\Delta G_{\mathit{vdW}}+\Delta G_{\mathit{hbond}}+\Delta G_{\mathit{desolv}}+\Delta G_{\mathit{elec}}+\Delta G_{\mathit{tor}}$$(1)

From the 1000 poses generated, only the binding pose with the most negative docking score was chosen for further analysis. The ligands and references were evaluated using the same scoring function. The docking scores were used to rank the compounds to determine which ones scored higher than the known inhibitors. The free energy of binding of each of the top poses of the 21 ligands was estimated using Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) methods. Docking experiments and the MM/PBSA calculation were done in triplicates.

Molecular dynamics simulations

Molecular dynamics (MD) simulations were carried out using Groningen Machine for Chemical Simulations (GROMACS) 2019.3 (Lemkul, 2018) on complexes with binding energies stronger than at least one of the references, estimated either by the docking scoring function or MMM/PBSA. The force field parameter utilized for the protein is the November 2018 version of Chemistry at Harvard Macromolecular Mechanics (CHARMM) 36 while the ligands’ force field parameters were generated from the CGenFF (https://cgenff.umaryland.edu/) server of the University of Maryland Baltimore (Vanommeslaeghe et al., 2010). A cubic box with periodic boundary conditions was used to contain the system, 10 Å from the protein’s edge. The box was then solvated using three-site transferable water model (TIP3P), to which Na⁺ ions were added to neutralize the system. Additional 0.15 M Na⁺ and Cl^- ions were added to mimic the physiological conditions. A 5,000-step structure optimization with steepest descent algorithm was performed to minimize the energy of the system, followed by a 100 ps equilibration in the NVT and NPT ensembles using V-rescale thermostat and Berendsen barostat for temperature and pressure coupling. The cutoff for both van der Waals and electrostatic interactions were kept at 1.2 nm, and the Particle Mesh Ewald method was used to calculate the long-range electrostatic interactions. In the production run, the temperature and pressure was kept at 310 K and 1 bar (Figures S2 and S3, Supplementary Materials). V-rescale thermostat with a time constant of 0.1 ps was used, while pressure was controlled by the Berendsen barostat using a time constant of 1 ps and a compressibility of 4.5 × 10⁻⁵bar⁻¹. Simulation run was 100 ns and snapshots of the trajectory were saved every 100 ps. All MD simulations were done in triplicates. The stability of the resulting trajectories were then evaluated using the root-mean-square deviation (RMSD), root-mean-square fluctuations (RMSF), hydrogen bond interactions, and MM/PBSA tools of GROMACS. In the H-bond analysis, the cutoff angle and distance were set to 30^° and 3.5 nm, respectively.

Figure 2. Reside interactions within 4.0 Å of the top eleven ligands and references with the SGLT-2 protein. References: (A) CAN, (B) DAP, (C) EMP. Candidate ligands: (D) AST, (E) BET, (F) CAS, (G) CEN, (H) ISO, (I) MAA, (J) MAD, (K) MEB, (L) MYR, (M) QUE, (N) RUT. Majority of the interactions were hydrophobic and polar while only a few electrostatic contacts were observed. Three residues (Glu 68, Ile 427, and Tyr 263) were also observed to have interacted with all the candidate ligands and the references.

Figure 3. Docking poses of the top five ligands and the references inside the binding pocket of SGLT-2. Also shown are the some of the interacting residues and their corresponding interaction types: hydrogen bonds (pink); aromatic hydrogen bonds (blue); pi-pi stacking (black). References: (A) CAN, (B) DAP, (C) EMP. Candidate ligands: (D) AST, (E) BET, (F) CAS, (G) CEN, (H) ISO, (I) MAA, (J) MAD, (K) MEB, (L) MYR, (M) QUE, (N) RUT. Extracellular (Phe 424, Tyr 87, Met 73) and intracellular (Tyr 263) gate residues are represented with green and purple spheres, respectively. Gln 428, the residue with the most critical role in the absorption mechanism, is represented with orange spheres.

Molecular mechanics/Poisson Boltzmann surface area

The Molecular Mechanics/Poisson Boltzmann method gives a more accurate estimation for the free energy of binding than the scoring function of docking experiments (Genheden & Ryde, 2015). It incorporates three energetic terms, including the molecular mechanics energies (bonded and non-bonded interactions), polar solvation energy, and non-polar solvation energy to solve for the change in the free energy. $$\Delta G=\Delta E_{\mathit{bonded}}+\Delta E_{\mathit{nonbonded}}+\Delta G_{\mathit{polar}}+\Delta G_{\mathit{nonpolar}}\mathrm{ }\left(2\right)$$

The g_mmpbsa tool (Kumari et al., 2014), a GROMACS compatible package, is used to implement the MM/PBSA calculations, and the python script MmPbSaStat.py provided in the g_mmpbsa package, compiles and estimates the interaction free energies. The tool currently does not include the calculation of the entropic energy and is therefore unable to give the absolute binding energy (Kumari et al., 2014).

The bonded interactions consisting of bond, angle, dihedral, and improper interactions of both the ligand and protein in the bound and unbound forms are assumed to not change so the first term of Equation (2) is zero. Non-bonded interactions include both electrostatic and van der Waals energies, which are based on the molecular mechanics force-field parameters. The last two terms, $\Delta G_{\mathit{polar}}$ and $\Delta G_{\mathit{nonpolar}}$ are the energies required to transfer a solute from a vacuum into the solvent. The former is estimated from the Poisson − Boltzmann (PB) equation, while latter makes use of solvent accessible surface area (SASA) model to estimate the repulsive and attractive forces associated with the cavity formation and van der Waals interactions (Kumari et al., 2014). Thus, Equation (2) can be transformed into: $${\Delta G=\Delta E_{\mathit{vdW}}+\Delta E}_{\mathit{elec}}+\Delta G_{\mathit{polar}}+\Delta G_{\mathit{SASA}}\mathrm{ }\left(3\right)$$

The clustering tool of GROMACS was also used to group together similar conformational states that are within 2.0 Å RMSD of each other in the stable phase of the trajectory. The resulting structures represent the most probable binding conformation of the ligand inside the protein.

Results and discussion

ADMET and drug-likeness of ligands

From the 45 properties and descriptors of QikProp, the three references scored favorably in 41 of these tests (Table 1). Methyl asiatate (MEA) obtained the highest consolidated QikProp score of 44, followed by bayogenin (BAY), madasiatic acid (MAA), and methyl brahmate (MEB), all with 43 points. Betulinic acid (BET), MEA, and MEB also had the most negative predicted aqueous solubility (QPlogS) that is within the recommended range, which correlates to high oral bioavailability. Moreover, the partition values (QPlogPC16, QPlogPoct, QPlogPw) of these compounds indicate high lipophilicity resulting to good passive diffusion and high gastrointestinal absorption (Shaikh & Siu, 2016). In addition, MEB had a ‘high’ HOA classification and scored 95.84% in PHOA, the highest among all phytochemicals and control compounds. Based on these prediction results, it is possible that MEB is a more potent drug than the three known SGLT-2 (PDB ID: 2XQ2) inhibitors.

Three compounds, asiaticoside (AST), madecassoside (MAD), and rutin (RUT), were rated poorly in QikProp, with an aggregate score of 23, 25, and 30 respectively. They also each incurred three violations of the Ro5 guidelines. These ligands got a ‘low’ classification in HOA and obtained a zero in PHOA. They also did not meet the criteria set on the allowable number of hydrogen bond donor, hydrogen bond acceptor, and molecular weight due to their relatively large structures. In accordance with Lipinski’s Rule of Five, these three compounds will have low chances of good permeability and solubility when taken orally. However, it also important to consider that Ro5 is an observation made only on known orally administered synthetic drugs and the rule does not include natural or semi-natural produced drugs. It has been determined that only 51% of all US FDA-approved drugs are compliant to Lipinski’s Rule, and macrocyclic and macrolide antibiotic drugs (i.e. erythromycin, cyclosporine, etc.), with molecular weight over 700 g/mol, obviously do not comply to these rules (Zhang & Wilkinson, 2007). But molecules with more than one violation should not be necessarily removed from further consideration but shall only receive low prioritization in the succeeding steps (Petit et al., 2012). RUT and AST have also been shown in in-vivo studies to possess anti-diabetic properties so it is worth studying whether the compounds’ glycemic control is a result of the inhibition of the SGLT-2 protein (Fitrianda et al., 2017; Supkamonseni et al., 2014). Moreover, among all of Centella asiatica’s phytochemicals, ASA, AST, MAD, and MAA are of primary interest since these are considered to be the plant’s most biologically active compounds (Hashim et al., 2011). Hence, even with poor oral bioavailability, AST, MAD and RUT were still included in the subsequent analysis. Other methods of drug administration could still be explored.

The online server, Molsoft (Molsoft, 2018) was also used to quantify the drug-likeness of a molecule. It calculates the drug-likeness score of the compounds using an algorithm based on the subject’s pharmacokinetic properties. These scores are then compared to the drug-likeness scores of known drugs (positives) and non-drugs (negatives) from their own set of databases. A drug-likeness score around 0 to 2 describes a drug-like molecule while a score between −6 to −1 corresponds to a non-drug compound. The drug-likeness scores of the majority of the ligands except for MYR are within the prescribed range (Table 1), indicating drug-like characteristics that are similar to most of the commercially approved drugs.

Molecular docking

The binding conformation and orientation of the phytochemicals were virtually the same in all three experiments for each ligand. From the 21 ligands, five exceeded the estimated binding energy of at least one of the references (Table 2). AST (–9.833 kcal/mol), MAD (–10.236 kcal/mol), and, RUT (–10.25 kcal/mol) had docking scores higher than all the reference drugs (–8.873 kcal/mol to −9.947 kcal/mol) suggesting that these three compounds have stronger binding affinities towards SGLT-2, while methyl brahmate (MEB) (–9.513 kcal/mol) and castilliferol (CAS) (–9.897 kcal/mol) only scored higher than one reference compound, DAP (–8.873 kcal/mol).

Table 2. Average docking and MM-PBSA scores of the 21 phytochemicals from the three trials.

Ligands		Molecular docking (kcal/mol)	MM-PBSA score (kcal/mol)
References
Canagliflozin	CAN	–9.9467 ± 0.0094	–12.1684 ± 0.1188
Dapagliflozin	DAP	–8.8733 ± 0.0946	–11.5946 ± 0.0471
Empagliflozin	EMP	–9.8367 ± 0.0929	–12.0542 ± 0.0154
Phytochemicals
Apigenin	API	–7.6600 ± 0.0356	–10.3889 ± 0.6390
Asiatic acid	ASA	–7.4233 ± 0.0377	–11.3328 ± 0.2318
Asiaticoside	AST	–9.9833 ± 0.0984	–13.4629 ± 0.2461
Bayogenin	BAY	–7.1267 ± 0.3159	–11.3511 ± 0.1928
Betulinic acid	BET	–8.5800 ± 0.0141	–14.6999 ± 0.4734
Brahmol	BRA	–8.0667 ± 0.3389	–9.1578 ± 0.2482
Castilliferol	CAS	–9.8967 ± 0.1173	–12.4843 ± 0.0434
Centellasapogenol	CEN	–7.8433 ± 0.2164	–13.3943 ± 0.2184
Isothankunic acid	ISO	–8.1600 ± 0.1344	–14.0346 ± 0.3757
Kaempferol	KAM	–7.9167 ± 0.1320	–9.7023 ± 0.7204
Madasiatic acid	MAA	–7.9167 ± 0.3795	–12.4086 ± 0.0807
Madecassic acid	MAC	–6.8600 ± 0.0927	–9.7708 ± 0.1017
Madecassoside	MAD	–10.2367 ± 0.1096	–14.3335 ± 1.0095
Methyl asiatate	MEA	–6.8567 ± 0.0047	–11.0141 ± 0.2290
Methyl brahmate	MEB	–9.5133 ± 0.0189	–12.8606 ± 0.2996
Myricetin	MYR	–8.2700 ± 0.0294	–12.4607 ± 0.2996
Patuletin	PAT	–7.2733 ± 0.1342	–8.4882 ± 0.6881
Pomolic acid	POM	–7.4767 ± 0.0471	–10.3767 ± 0.5551
Quercetin	QUE	–8.5867 ± 0.1281	–13.5341 ± 0.8829
Rutin	RUT	–10.2500 ± 0.1395	–15.3976 ± 0.5458
Terminolic acid	TER	–6.8533 ± 0.9595	–8.6702 ± 0.6397

Eleven compounds had stronger MM/PBSA binding free energies than at least one of the references. The top five compounds with the highest docking scores still obtained stronger MM/PBSA binding energies than the references (–11.595 kcal/mol to −12.168 kcal/mol). RUT also still obtained the highest energy (–15.266 kcal/mol) from all phytochemicals tested. We can observe from the results that the MM/PBSA energies do not necessarily follow the docking scores. For example, BET got the second highest MM/PBSA energy (–14.700 kcal/mol) but was not part of the top compounds with the best docking scores. This is also apparent on some compounds including centellasapogenol (CEN), isothankunic acid (ISO), madasiatic acid (MAA), myricetin (MYR), and quercetin (QUE), which all had favorable MM/PBSA energies. The difference in the scores could be due to the inclusion of polar and non-polar solvation energies which resulted to a more negative binding energy in the MM/PBSA method (Thompson et al., 2008). Only the ligands with a stronger binding affinity than at least one of the references, evaluated either by the docking scoring function or MM/PBSA, were included in the succeeding experiments. These are AST, BET, CAS, CEN, ISO, MAA, MAD, MEB, MYR, QUE, and RUT.

Ligand – receptor interactions

The majority of the residues in the binding site of SGLT-2 that interacted with the three references, also interacted with the eleven candidate ligands (Table 3). Most of interactions are hydrophobic and polar in nature, with only a few electrostatic interactions were present (Figure 2). Hydrophobic interactions are the most common type of interactions in protein-ligand complexes, and usually determine the overall stability of a drug-receptor interactions (de Freitas & Schapira, 2017). Gln 428 and the external gate residue Tyr 263, which both hold critical roles in binding activity, were observed to form polar and hydrophobic interactions with the majority ligands, respectively (Figure 2 and Table 3). In contrast, there were no visible interaction with the gate residue Met 73 in any of the compounds, but interactions with the other extracellular gate residues, such as Tyr 87 and Phe 424, were present in some (Figure 3). MAD was the only ligand to miss a contact with any of the external gate residues. Electrostatic interaction with Glu 68 was also observed to be consistent with all the ligands and the references. Aside from its electrostatic contribution, Glu 68 also produced hydrogen bonding (H-bonding) interactions with all the eleven candidate ligands (Table 4). H-bonding have been shown in studies to have significant importance in conferring stability to protein-ligand interactions (Li et al., 2018; Wang et al., 2011). Glu 68 produced three H-bonds with ISO, two H-bonds with CAN, DAP, MYR, QUE, and RUT, and one H-bond with the rest of the ligands except for AST and MEB – the only two compounds that Glu 68 failed to pair with (Table 4). H-bond contacts were also present on other residues including Ala 63, Asn 260, Asn 267, Asp 189, Gln 428, Glu 88, Ser 66, Ser 368, Thr 431, Tyr 138, while CEN, MYR, and QUE formed H-bond connections with Gln 428 – one of the important residues.

Table 3. Common interacting residues among the ligands with the highest docking scores and the references within 4.0 Å.

Residues	References				Candidate ligands
	CAN	DAP	EMP	AST		BET	CAS	CEN	ISO	MAA	MAD	MEB	MYR	QUE	RUT
Ala 63	✓	✓	✓	✓			✓	✓	✓	✓		✓	✓	✓	✓
Ala 259	✓	✓	✓										✓	✓	✓
Asn 64	✓	✓	✓				✓	✓				✓	✓	✓	✓
Asn 142	✓			✓		✓	✓		✓	✓	✓	✓			✓
Asn 260	✓	✓	✓	✓		✓						✓	✓	✓	✓
Asn 267	✓			✓					✓	✓	✓				✓
Gln 69	✓	✓	✓			✓	✓	✓		✓	✓		✓	✓	✓
Gln 268											✓				✓
Gln 428	✓	✓	✓	✓		✓	✓	✓		✓	✓	✓	✓	✓	✓
Glu 68	✓	✓	✓	✓		✓	✓	✓	✓	✓	✓	✓	✓	✓	✓
Glu 88	✓						✓	✓					✓	✓	✓
Ile 427	✓	✓	✓	✓		✓	✓	✓	✓	✓	✓	✓	✓	✓	✓
Met 369			✓	✓		✓		✓	✓	✓		✓
Phe 424		✓	✓			✓	✓	✓		✓		✓	✓	✓	✓
Ser 66	✓	✓	✓	✓			✓	✓	✓				✓	✓	✓
Ser 91		✓		✓			✓					✓	✓	✓	✓
Ser 368				✓				✓	✓	✓
Ser 435		✓				✓
Thr 431	✓	✓	✓				✓		✓		✓				✓
Trp 264		✓	✓				✓								✓
Tyr 87	✓	✓	✓	✓		✓	✓					✓	✓	✓	✓
Tyr 138				✓				✓	✓	✓					✓
Tyr 262	✓	✓	✓			✓			✓		✓
Tyr 263	✓	✓	✓	✓		✓	✓	✓	✓	✓	✓	✓	✓	✓	✓
Trp 264		✓	✓			✓	✓						✓	✓	✓

Three residues, namely Glu 68, Ile 427, and Tyr 263, were identified to have interacted with all the ligands. Among these, Glu 68 and Tyr 263 are known to play critical roles in the binding activity of the protein.

Table 4. H-bond interactions of the top candidate ligands.

Ligand	No. of H-bonds	Interacting residue (H-bond length, Å)
References
CAN	4	Glu 68 (1.65), Glu 68 (1.92), Ala 63 (2,10), Tyr 263 (2.15)
DAP	3	Glu 68 (1.82), Glu 68 (2.01), Ser 66 (2.00)
EMP	3	Glu 68 (1.87), Ser 66 (2.50), Ala 63 (1.97)
Phytochemicals
AST	2	Asn 260 (2.48), Asp 189 (2.46)
BET	1	Asn 260 (2.35)
CAS	3	Glu 68 (1.75), Glu 88 (2.05), Asn 142 (2.09)
CEN	2	Gln 428 (2.46), Ala 63(1.96)
ISO	4	Glu 68 (2.14), Glu 68 (2.17), Glu 68 (2.15), Thr 431 (1.96)
MAA	3	Gln 428 (2.01), Ser 368 (2.44), Tyr 138 (1.85)
MAD	2	Glu 68 (2.31), Asn 267 (2.38)
MEB	1	Ala 63 (1.89)
MYR	4	Asn 260 (1.82), Gln 428 (1.99), Glu 68 (2.03), Glu 68 (1.80)
QUE	4	Asn 260 (1.79), Gln 428 (1.95), Glu 68 (2.12), Glu 68 (1.84)
RUT	4	Asn 260 (1.87), Glu 68 (2.02), Glu 68 (2.14), Thr 431 (1.77)

Of particular interest is the residue Glu 68. It produced H-bonds connections with almost all ligands and had multiple contacts with CAN, DAP, ISO, MYR, QUE, and RUT.

Molecular dynamics simulations

To evaluate the stability of the resulting complexes from molecular docking, a time dependent MD simulation was conducted for 100 ns using GROMACS 2019.3. Two types of RMSD were used to quantify how much the complex has changed over the simulation time. Complex-RMSD (Figure 4) describes the movement of both the ligand and the protein backbone while ligand-RMSD (Figure 5) shows the fluctuation of the ligand alone inside the binding pocket – both of which are used to describe the stability. Table 5 tabulates the average complex-RMSD and ligand-RMSD on the last 10 ns of the simulation.

Figure 4. Complex-RMSD of the top five ligands and the three references docked against the binding site of SGLT-2. References: (A) CAN, (B) DAP, (C) EMP. Candidate ligands: (D) AST, (E) BET, (F) CAS, (G) CEN, (H) ISO, (I) MAA, (J) MAD, (K) MEB, (L) MYR, (M) QUE, (N) RUT.

Figure 5. Ligand-RMSD of the top five ligands and the three references docked against the binding site of SGLT-2. References: (A) CAN, (B) DAP, (C) EMP. Candidate ligands: (D) AST, (E) BET, (F) CAS, (G) CEN, (H) ISO, (I) MAA, (J) MAD, (K) MEB, (L) MYR, (M) QUE, (N) RUT.

Table 5. Complex-RMSD and Ligand-RMSD of the top eleven promising compounds.

	Complex-RMSD				Ligand-RMSD
	Min	Max	Mean	Standard deviation	Min	Max	Mean	Standard deviation
References
CAN	2.8416	3.1596	3.0224	0.1334	1.2726	1.8252	1.5685	0.2273
DAP	2.7838	3.1376	2.9177	0.1567	1.4361	2.1042	1.8371	0.2887
EMP	2.5993	2.6496	2.6232	0.0206	1.8895	2.0450	1.9634	0.0637
Phytochemicals
AST	3.1188	3.3321	3.2164	0.0880	2.4405	2.7547	2.6384	0.1407
BET	2.4166	3.3621	2.8874	0.2263	0.3464	0.7171	0.5258	0.0351
CAS	2.5666	2.8491	2.6689	0.1278	1.0802	1.8342	1.5278	0.3236
CEN	2.2290	3.0819	2.6854	0.1674	1.0366	1.4649	1.2570	0.0394
ISO	2.2290	3.1154	2.6186	0.1990	0.3088	0.9537	0.6204	0.1580
MAA	2.2747	2.8602	2.5220	0.0187	1.0174	1.4532	1.2426	0.0455
MAD	3.2872	3.3551	3.3158	0.0287	2.4349	2.6266	2.5605	0.0888
MEB	2.3618	2.7272	2.5007	0.1616	1.3354	1.4208	1.3708	0.0363
MYR	2.9423	3.7539	3.3814	0.1134	0.5220	1.5905	1.3091	0.0960
QUE	3.1137	3.8123	3.4176	0.0093	0.2285	1.3705	0.6417	0.1947
RUT	2.4970	2.6477	2.5498	0.0693	1.4426	2.2215	1.9320	0.3480

Also shown are the RMSDs of the three references, CAN, DAP and EMP. Tabulated values are obtained from the last 10 ns of the simulation.

All complexes except for AST and MAD equilibrated within the simulation time of 100 ns, achieving an RMSD range between 1 to 3 Å. CAS, CEN, MAA, MEB, and RUT stabilized the quickest at around 40 ns (Figure 4), while BET, MYR, and QUE equilibrated at around 60 ns. ISO took a longer time before it reached equilibrium at 80 ns in two of the three runs.

If we look at the average number of H-bond pairs from all three runs (Figure 6), CEN, MYR, MEB and QUE, had consistent number of H-bonds throughout the simulation. MEB even had a gradual increase of H-bonding towards the end which further conferred the rigidity of the complex. BET, CAS, and MAA may have fewer H-bonding, but these molecules are smaller, so the number of its H-bond contacts were just enough to hold the ligands in place. ISO formed the least number of H-bond pairs that may have affected the complex’s stabilization. Meanwhile, RUT’s H-bonds connections sharply decreased at the 40 ns mark but remained consistent afterwards.

Figure 6. Average number of intermolecular hydrogen bond interactions from the three runs throughout the entire simulation. References: (A) CAN, (B) DAP, (C) EMP. Candidate ligands: (D) AST, (E) BET, (F) CAS, (G) CEN, (H) ISO, (I) MAA, (J) MAD, (K) MEB, (L) MYR, (M) QUE, (N) RUT.

AST and MAD had increasing complex-RMSD (Figure 4) towards the latter part of the simulation, and large and erratic fluctuations were also observed in their ligand-RMSD (Figure 5). The number of H-bonds contacts of AST continued to slope downwards (Figure 6), and while MAD registered a steady number of H-bonds pairs, it was relatively low considering its size so it may have needed more H-bond connections to fully secure itself inside the binding pocket of SGLT-2. The large structures of AST and MAD may have caused the ligands to form interactions outside the binding region.

RMSF plots are then used to describe the flexibility of each amino acid residue during molecular dynamics simulations (Benson & Daggett, 2012). It was observed that the RMSF values of the residues constituting the active site of SGLT-2 were only slightly changed (Figure 7). The low values at the troughs describe the residues’ reduced movements which prevented the binding site to undergo drastic structural changes and made it less susceptible to rearrangement. Since the shape of the site was not altered during the simulation, this allowed the ligands to stay and fully latch themselves inside the cavity. The brown bars on the x-axis represent the position of the critical residues including Gln 428, and the gate residues, Met 73, Tyr 263, Phe 424, and Tyr 87. As seen in Figure 7, all these residues had low RMSF values. It follows that the ligands that sufficiently formed contacts with the aforementioned active site residues also had constrained movements.

Figure 7. RMSF plot of the top five ligands and the three references docked against the binding site of SGLT-2. References: (A) CAN, (B) DAP, (C) EMP. Candidate ligands: (D) AST, (E) BET, (F) CAS, (G) CEN, (H) ISO, (I) MAA, (J) MAD, (K) MEB, (L) MYR, (M) QUE, (N) RUT. The brown bars in the x axis represent the location of critical residues including Gln 428, and the gate residues, Met 73, Tyr 263, Phe 424, and Tyr 87.

MM/PBSA analysis

The binding affinities (Equation (3)) of the eleven complexes were estimated from the average MM/PBSA binding free energies using the last 10 ns (last 100 frames) of the trajectories.

Five ligands (AST, BET, CEN, MEB, RUT) obtained more negative MM/PBSA energies than at least one of the references. MEB (–31.07 kcal/mol) and BET (–29.60 kcal/mol) registered the strongest binding energies and surpassed the energies of all the references (EMP: –18.113 kcal/mol, DAP: –21.16 kcal/mol, CAN: –28.763 kcal/mol) (Table 6). CEN got the third highest binding free energy (–24.64 kcal/mol) but it was not able to match the energy of CAN (–28.763 kcal/mol). AST (–19.10 kcal/mol) and RUT (–18.59 kcal/mol) only managed to score higher than the reference EMP (–18.11 kcal/mol), while CAS, ISO, MAA, MAD, MYR, and QUE failed to surpass the binding energies of any of the reference compounds.

Table 6. Average binding free energies of the top ligands.

Ligand	van der Waals energy (kcal/mol)	Electrostatic energy (kcal/mol)	Polar solvation energy (kcal/mol)	SASA energy (kcal/mol)	Binding energy (kcal/mol)
CAN	–49.69 ± 2.62	–9.27 ± 2.80	35.86 ± 3.97	–5.75 ± 0.21	–28.76 ± 3.47
DAP	–42.49 ± 3.08	–8.73 ± 4.16	35.37 ± 4.14	–5.31 ± 0.21	–21.16 ± 3.41
EMP	–41.50 ± 3.05	–13.16 ± 3.37	42.08 ± 4.41	–5.53 ± 0.22	–18.11 ± 3.81
AST	–72.27 ± 4.12	–29.50 ± 5.95	92.09 ± 12.27	–9.42 ± 0.41	–19.10 ± 9.88
BET	–49.64 ± 2.62	–4.34 ± 1.95	29.85 ± 3.60	–5.47 ± 0.21	–29.60 ± 3.86
CAS	–43.33 ± 2.22	–24.48 ± 3.95	54.89 ± 4.20	–5.54 ± 0.23	–17.44 ± 3.39
CEN	–42.98 ± 2.73	–7.84 ± 2.26	31.49 ± 2.87	–5.32 ± 0.23	–24.64 ± 2.60
ISO	–40.25 ± 3.31	–6.47 ± 3.16	34.91 ± 4.74	–5.14 ± 0.24	–16.95 ± 3.72
MAA	–40.28 ± 2.34	–5.41 ± 2.20	33.42 ± 3.72	–5.21 ± 0.20	–17.48 ± 3.41
MAD	–75.99 ± 4.56	–36.34 ± 6.51	104.35 ± 7.90	–9.76 ± 0.31	–17.73 ± 5.56
MEB	–52.51 ± 2.56	–8.75 ± 2.02	36.10 ± 2.88	–5.92 ± 0.21	–31.07 ± 3.32
MYR	–31.20 ± 3.52	–35.86 ± 3.85	58.61 ± 3.33	–4.03 ± 0.15	–12.47 ± 3.33
QUE	–30.66 ± 2.95	–25.88 ± 3.72	43.42 ± 4.14	–3.81 ± 0.17	–16.93 ± 3.58
RUT	–49.10 ± 3.28	–35.61 ± 4.38	72.66 ± 4.76	–6.54 ± 0.22	–18.59 ± 4.08

The binding energy is estimated from the sum of the van der Waal, electrostatic, polar solvation, and solvent accessible surface area (SASA) energies during the last 10 ns of the simulation.

The magnitude of the energetic components of MEB closely resembles the values obtained by the two references, CAN and DAP. The strong van der Waals interactions of MEB, BET, and CEN compensated for their weak electrostatic forces that led to a considerably high binding affinity. ISO and MAA had similar electrostatic, polar solvation, and SASA energies as MEB, BET, and CEN, but their van der Waals energies were not as strong. It was also observed that molecules with higher polar surface area (Table S2 (PSA), Supplementary Materials) such as AST, MYR, and QUE, and RUT tend to have higher polar solvation energies but more negative electrostatic energies. Polar solvation energy, along with the non-polar solvation energy, constitutes the energy required to transfer a solute from a vacuum into the solvent. It describes the work done for moving a partial charge of a solute in a reference medium to another point in a solvent medium (Izairi & Kamberaj, 2017; Kumari et al., 2014). It is expected that molecules with large dipole moments or high polar surface area (PSA) would require more work in the desolvation process. This could explain the high polar solvation energies of AST, MAD, MYR, QUE, and RUT. BET, CEN, and MEB are not as polar as all the other ligands (Table S2 (PSA), Supplementary Materials) and only required relatively low polar solvation energies in the coupling process (Table 6). RUT obtained one of the most negative electrostatic energies, which can be correlated to its high dipole moment – the highest among all ligands. When all the energetic terms are totaled, RUT was still able to surpass the binding energy of one of the references, EMP, hence RUT can also be considered as a potential candidate molecule.

Decomposition of the binding free energy on a per residue basis of the ligands with the strongest MM/PBSA energies (AST, BET, CEB, MEB, RUT) is shown in Figure 8, while the decomposition of top contributing residues into energy contributions from molecular mechanics energy, polar solvation energy, and non-polar solvation energy is presented in Figure 9. The residues Tyr 263, Phe 424, Met 369, and Ile 427 were among the highest contributors to the observed negative binding free energies (Figure 8). It is also notable to add that Tyr 263 and Phe 424 are part of the gate residues, which are important in the inhibition mechanism of SGLT-2. Tyr 263 consistently imparted the strongest molecular mechanics energies which largely compensated the unfavorable positive solvation energy (Figure 9). Residues Tyr 263, Phe 424, Met 369, and Ile 427 also contributed positive polar solvation energies which favors complex dissociation. Nevertheless, the total energy contribution from the three energetic terms still favored complex formation driven by the strong intermolecular interactions between the ligand and the protein. On the other hand, Asp 336 and Gln 429 provided the most negative polar solvation energies in the RUT-receptor complex. Asp 336 and Gln 429 also imparted favorable polar solvation energies in the MEB-receptor and BET-receptor complexes, respectively.

Figure 8. Decomposition of the binding free energy on a per residue basis of the complexes with the strongest MM/PBSA energies. References: (A) CAN, (B) DAP, (C) EMP. Final candidate ligands: (D) AST, (E) BET, (F) CEN, (G) MEB, (H) RUT.

Figure 9. Residues with the highest energy contribution from the complexes with the strongest MM/PBSA energies. Also shown is the decomposition of the binding free energy on a per-residue basis into contributions from the molecular mechanics (van der Waals + electrostatic) energy, polar solvation energy, and non-polar solvation energy. References: (A) CAN, (B) DAP, (C) EMP. Final candidate ligands: (D) AST, (E) BET, (F) CEN, (G) MEB, (H) RUT.

Geometric clustering within 2.0 Å RMSD was done on the last 10 ns of the simulations to generate representative structures of each complex (Figure 10). These representative structures were used to determine ligand-residue contacts. BET, CEN, ISO, MEB, QUE, retained the hydrophobic and polar interactions with the important residues that were also present in their initial docking pose (Figure 3). BET and QUE kept their interactions with Phe 424, Tyr 87, Tyr 263, and Gln 428; CEN and MEB with Phe 424, Tyr 263, and Gln 428; and ISO with Tyr 263. AST and RUT were unsuccessful to sustain an interaction with Tyr 87; nonetheless, AST formed new interactions with two external gate residues, Met 73 and Phe 424. MEB also established a hydrophobic interaction with Met 73, the external gate residue that did not form any connections with any of the ligands before the MD simulations. MYR lost its interaction with Gln 428 but still interacted with three gate residues. CAS, MAA and MAD also lost their initial contacts with all the important residues except Tyr 263, but MAD made a new connection with Phe 424. Interactions with important residues Tyr 263, Phe 424, and Gln 428 were also noticeably present in the ligands with the strongest MM/PBSA energies (AST, BET, CEN, MEB, and RUT). These observations point towards the high potential of AST, BET, CEN, MEB, and RUT as inhibitors of SGLT-2.

Figure 10. Representative structures of the most populated cluster of conformations that are within 2.0 Å RMS of each other in the last 10 ns of the trajectory: hydrogen bonds (pink); aromatic hydrogen bonds (blue); pi-pi stacking (black). References: (A) CAN, (B) DAP, (C) EMP. Candidate ligands: (D) AST, (E) BET, (F) CAS, (G) CEN, (H) ISO, (I) MAA, (J) MAD, (K) MEB, (L) MYR, (M) QUE, (N) RUT. Extracellular (Phe 424, Tyr 87, Met 73) and intracellular (Tyr 263) gate residues are represented with green and purple spheres, respectively. Gln 428, the residue with the most critical role in the absorption mechanism, is represented with orange spheres.

Conclusion

Molecular docking experiments were performed on 21 phytochemicals from Centella asiatica, a known folk-medicinal plant, and three reference drugs, onto the SGLT-2 receptor. This revealed 11 compounds (AST, BET, CAS, CEN, ISO, MAA, MAD, MEB, MYR, QUE, and RUT) with better docking scores or MM/PBSA energies than at least one of the references. 100-ns MD simulations were done on these 11 complexes to determine their stability, binding energies, and interacting residues. RMSD and RMSF trajectory analyses elucidated the stable behavior of the ligands in the binding region. MM/PBSA free energy calculations using the last 10 ns of the simulations showed that AST, BET, CEN, MEB, and RUT had stronger binding free energies than at least one of the references. Compounds with higher polar surface area tend to have weaker affinity towards the protein because of the larger polar solvation energy requirement during coupling. Detailed interaction analysis revealed the importance of the residues Tyr 263, Phe 424, Met 369, and Ile 427 for their molecular mechanics energy contributions, and Asp 336 and Glu 429, for their favorable polar solvation energies.

All these demonstrate the potential of AST, BET, CEN, MEB, and RUT as SGLT-2 inhibitors. The five candidates matched the binding affinities of known SGLT-2 inhibitor drugs, and formed favorable interactions with known critical residues of the receptor. Among the five, however, BET, CEN, and MEB have better absorption, distribution, metabolism, excretion properties and better equilibration profiles than AST and RUT. Their trajectories inside SGLT-2 are also comparable to those of the reference molecules. Thus, further in vitro and in vivo studies using BET, CEN and MEB as inhibitors against the SGLT-2 protein are recommended.

Abbreviations
ADMET	=	absorption, distribution, metabolism, excretion, toxicity
ASA	=	asiatic acid
AST	=	asiaticoside
BET	=	betulinic acid
BAY	=	bayogenin
CAN	=	canagliflozin
CEN	=	centellasapogenol
CAS	=	castilliferol
DAP	=	dapagliflozin
EMP	=	empagliflozin
GROMACS	=	Groningen Machine for Chemical imulations
ISO	=	isothankunic acid
MAA	=	madasiatic acid
MAD	=	madecassoside
MD	=	molecular dynamics
MEA	=	methyl asiatate
MEB	=	methyl brahmate
MM/PBSA	=	molecular mechanics Poisson–Boltzmann surface area
MYR	=	myricetin
Ro5	=	Lipinski’s Rule of Five
QUE	=	quercetin
RUT	=	rutin
SGLT-2	=	Sodium-glucose co-transporter 2

Acknowledgments

The authors want to thank the Computing and Archiving Research Environment of the Department of Science and Technology – Advanced Science and Technology Institute (COARE-DOST-ASTI), the Engineering Research and Development for Technology of the Department of Science and Technology (ERDT-DOST), and the University of the Philippines Diliman – Office of the Vice-Chancellor for Research and Development (UP-OVCRD) for their support to this research.

Disclosure statement

The authors have no conflicts of interest to report.