*Article* **Rational Design of Novel Inhibitors of** α**-Glucosidase: An Application of Quantitative Structure Activity Relationship and Structure-Based Virtual Screening**

**Sobia Ahsan Halim <sup>1</sup> , Sumaira Jabeen 2 , Ajmal Khan 1, \* and Ahmed Al-Harrasi 1, \***


**Abstract:** α-Glucosidase is considered a prime drug target for Diabetes Mellitus and its inhibitors are used to delay carbohydrate digestion for the treatment of diabetes mellitus. With the aim to design α-glucosidase inhibitors with novel chemical scaffolds, three folds ligand and structure based virtual screening was applied. Initially linear quantitative structure activity relationship (QSAR) model was developed by a molecular operating environment (MOE) using a training set of thirty-two known inhibitors, which showed good correlation coefficient (r <sup>2</sup> = 0.88), low root mean square error (RMSE = 0.23), and cross-validated correlation coefficient r 2 (q <sup>2</sup> = 0.71 and RMSE = 0.31). The model was validated by predicting the biological activities of the test set which depicted r <sup>2</sup> value of 0.82, indicating the robustness of the model. For virtual screening, compounds were retrieved from zinc is not commercial (ZINC) database and screened by molecular docking. The best docked compounds were chosen to assess their pharmacokinetic behavior. Later, the α-glucosidase inhibitory potential of the selected compounds was predicted by their mode of binding interactions. The predicted pharmacokinetic profile, docking scores and protein-ligand interactions revealed that eight compounds preferentially target the catalytic site of α-glucosidase thus exhibit potential αglucosidase inhibition *in silico*. The α-glucosidase inhibitory activities of those Hits were predicted by QSAR model, which reflect good inhibitory activities of these compounds. These results serve as a guidelines for the rational drug design and development of potential novel anti-diabetic agents.

**Keywords:** α-Glucosidase; QSAR modeling; homology modeling; molecular docking; ADMET profiling

#### **1. Introduction**

Diabetes mellitus (DM) is one of the most prominent metabolic, chronic and persistent diseases illustrated by hyperglycemia, which itself poses a great global health challenge [1,2]. The morbidity rate of DM is expected to exceed up to 10.4% by 2040 globally [3]. DM is linked with impaired production of insulin or due to improper response of body cells towards insulin. Several complications including hypertension, kidney failure, loss of vision, neuropathy, cardiovascular diseases (CVDs) and atherosclerosis are associated with inadequate production of insulin or due to prolonged hyperglycemia. DM is classified into Type I DM (T1DM) and Type II DM (T2DM). T1DM is also known as insulin dependent DM, which is caused by environmental, genetic and autoimmune factors. T-cell mediated autoimmune damage of *β*-cells of the pancreas is responsible for this type and illustrated by complete scarcity of insulin. One of the significant complications of this disease is ketoacidosis. Patients are susceptible to other diseases like Grave's disease, Addison's disease, Celiac sprain and autoimmune hepatitis [3]. T2DM is non-insulin dependent diabetes or adult onset diabetes associated with insulin resistance in peripheral tissues

**Citation:** Halim, S.A.; Jabeen, S.; Khan, A.; Al-Harrasi, A. Rational Design of Novel Inhibitors of α-Glucosidase: An Application of Quantitative Structure Activity Relationship and Structure-Based Virtual Screening. *Pharmaceuticals* **2021**, *14*, 482. https://doi.org/ 10.3390/ph14050482

Academic Editor: Osvaldo Andrade Santos-Filho

Received: 3 March 2021 Accepted: 2 April 2021 Published: 19 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

or relative deficiency of insulin. It is a progressive and most common type of disease. *β*-Cells are not damaged in this case; however, obesity is associated in most individuals and ketoacidosis rarely occurs. Most of the time, patients remain asymptomatic and have an enhanced risk of serious macrovascular and microvascular complications. Moreover, autooxidation of glucose and enhanced oxidative stress are also associated with T2DM. The risk of T2DM increases due to increased age, weight, improper physical activity and sedentary lifestyle [4]. In the advance stage of T2DM more complications occur as free radicals are formed, which cause insulin resistance and the dysfunction of *β*-cells.

Enhanced production of glucose is directed by α-glucosidase in the body; therefore, its level can be controlled by discovering its inhibitors [5,6]. α-Glucosidase is a lysosomal exo-glycosidase that catalyzes the breakdown of complex sugar like starch and disaccharides to glucose that are further absorbed by gut and subsequently raises postprandial hyperglycemia. Thus, α-glucosidase is the pathological hallmark of this disease. Furthermore, this enzyme has been linked in tumor metastasis in association with collagen type I and IV. In DM patients, inhibition of α-glucosidase efficiently lowers the risk of colorectal cancer and cerebrovascular events [7,8].

Many drugs like sulfonylureas, meglitinides, biguanides and thiazolidinediones and α-glucosidase inhibitors are most potent hypoglycemic agents [9]. Metformin, miglitol, acarbose, voglibose and glibenclamide are clinically used drugs for the management of diabetes. They are highly expensive and put a financial burden on patients, which results in huge economic loss. Their detrimental side effects also created an alarming situation like dropsy, increase in weight, resistance against drug, liver disorder, renal tumors, acute hepatitis, hepatic injury and hypoglycemia [10]. Thus, design and discovery of novel α-glucosidase inhibitors with reduced cost and lower side effects are urgently needed to control DM [11–13].

Computational drug designing approaches have proven effective in delivering novel chemical scaffolds in the market against several drug targets [14–16]. Due to our interest in finding novel leads against α-glucosidase [17–19], ligand and structure based virtual screening strategies were applied in this study to target α-glucosidase enzyme. The study comprises of ligand based QSAR modeling, structure-based virtual screening, pharmacokinetic profiling, and QSAR-based prediction of biological activities of new scaffolds against α-glucosidase. In-silico techniques have proven effective in the discovery of chemical agents against specific drug target with high binding affinity and specificity [15]. Several potent drugs including Captopril (Angiotensin-converting enzyme inhibitor, treats congestive heart failure and hypertension), Nelfinavir (HIV-protease inhibitor), Zanamivir (neuraminidase inhibitor of influenza viruses), Dorzolamide (carbonic anhydrase inhibitor), Cinanserin (potent inhibitor of 3C-like protease of SARS), Saquinavir, Indinavir, and Ritonavir (anti-HIV/AIDS drugs) were successfully designed by computational methods which are available on the market [20–23]. In addition, computer-aided drug designing (CADD) has delivered numerous novel anti-cancer compounds such as Erlotinib (EGFR kinase inhibitor, suppresses Non-small cell lung carcinoma and pancreatic cancer), Sorafenib (VEGFR inhibitor, treat thyroid, renal and liver cancer), Lapatinib (for the management of ERBB2-postive breast cancer), Abiraterone (inhibits androgen production to treat metastatic castration-resistant prostate cancer or hormone-refractory prostate cancer) and Crizotinib (ALK inhibitor, to treat Non-small cell lung carcinoma) [24–39].

Keeping the successful applications of CADD in mind, we employed both ligand and structure-based methods to design novel α-glucosidase inhibitors.

#### **2. Results and Discussion**

#### *2.1. QSAR Modeling*

QSAR models are widely used to predict the biological activities of compounds in silico [40,41]. Molecular Operating Environment (MOE) use descriptor-based technique to develop linear QSAR model and give two values, i.e., predicted pIC<sup>50</sup> and residual value. Initially 192 2D-descriptors were calculated for thirty-two known α-glucosidase inhibitors

(α**GI1-**α**GI4,** α**GI6-**α**GI15,** α**GI17-**α**GI20,** α**GI22-**α**GI24,** α**GI26-**α**GI31,** α**GI33-**α**GI34,** α**GI36-**α**GI38**) via QuaSAR-Descriptor module and QuaSAR-contingency was used to select the appropriate descriptors. Out of 192, nine descriptors (apol, bpol, a\_acc, a\_heavy, logP(o/w), logS, TPSA, Weight and wienerPol) were suggested as best descriptors for our data by contingency analysis. Partial least square (PLS) regression method was applied on training set to develop model. Model adequacy was measured as the square of correlation coefficient (r<sup>2</sup> ), root mean square error (RMSE), cross–validated r<sup>2</sup> (r2 LOO or q<sup>2</sup> ) and cross– validated RMSE. The developed model showed good correlation coefficient (r<sup>2</sup> = 0.88) and low root mean square error (RMSE = 0.23), suggesting the strength of the model. Cross validation was performed using leave-one-out (LOO) cross validations scheme which resulted acceptable cross-validated correlation coefficient (q<sup>2</sup> = 0.71 and RMSE = 0.31). The biological activities of the test set [six molecules (α**GI5,** α**GI16,** α**GI21,** α**GI25,** α**GI32,** α**GI35**)] were predicted by the model which showed r<sup>2</sup> value of 0.82. The model accurately predicted the activities of test set. The chemical structures, biological activities, predicted activities and residual values of α**GI**s are summarized in Table 1. The correlation plots of training and test set are shown in Figure 1. The 2D linear regression model is shown in the equation below, which is statistically significant and explain 88% of the variability of the IC<sup>50</sup> coefficient, characterized by usefulness of the model to predict the α-glucosidase inhibitory activity of test set compounds.

**Table 1.** The chemical structures, α-Glucosidase inhibitory activities, QSAR-predicted activities, and residual values of [α-Glucosidase known Inhibitors (α**GIs**)] used in QSAR modeling.



**Table 1.** *Cont*.


**Table 1.** *Cont*.


**Table 1.** *Cont*.

\* Test Set Compounds.

**Figure 1.** Graphical presentation of correlation between actual and predicted α-glucosidase inhibitory activities of training and test set compounds. R<sup>2</sup> of training and test set are 0.883 and 0.820, respectively.

One metric to measure a descriptor's relative importance in a QSAR model is to compute the product of the magnitude of a regression coefficient times the range in values adopted by its descriptor across the training set. This metric is a measure of the variability of dependent variable of the QSAR upon the descriptor. The relative importance of descriptors is: a\_heavy = 1.000, Weight = 0.773, and wienerPol = 0.718, TPSA = 0.471, a\_acc = 0.424, logP(*o/w*) = 0.355, logS = 0.244, apol = 0.242, bpol = 0.207. The data shows that the weight, and wienerPol descriptors are the most important factor because of their correlation coefficient > 0.5, followed by TPSA, a\_acc and logP.

The positive coefficient of descriptor in model reflects that increasing the number of heavy (polar) atoms, hydrogen bond acceptor atoms, logP(o/w), bpol (Sum of the absolute value of the difference between atomic polarizabilities of all bonded atoms in the molecule (including implicit hydrogens) with polarizabilities) and topological polar surface area (Å<sup>2</sup> ) will increase the activity of compounds. The descriptors a\_heavy, logP(o/w), bpol and TPSA are physical properties which can be calculated from the connection table (with no dependence on conformation) of a molecule while a\_acc belongs to pharmacophore feature/atom type descriptors.

The negative coefficient of descriptors shows that decreasing some of the physical properties including apol, logS, molecular weight and one of the adjacency and distance matrix descriptors, i.e., wiener polarity number will be beneficial to enhance the activity of the compounds.

$$\begin{array}{l} \text{pIC50} = 3.39574 \ -0.08965 \ \text{(apol)} \ + 0.09305 \ \text{(bpol)} \ + 0.49740 \ \text{(a.acc.)}\\ + 0.71639 \ \text{(a\\_heavy)} \ + 0.40979 \ \text{(logP(o/w))} \ - 0.37591 \ \text{(logS)}\\ + 0.02502 \ \text{(TPSA)} \ - 0.03659 \ \text{(Weight)} \ - 0.24818 \ \text{(wienerPol)} \end{array} \tag{1}$$

$$\text{r2} = 0.88, \text{ RMSE} = 0.23, \text{ q2} = 0.71, \text{ RMSE}(\text{LCO}) = 0.31, \text{ N(train)} = 32; \text{ N(test)} = 64$$

where apol = Sum of the atomic polarizabilities (including implicit hydrogens) with polarizabilities, bpol = Sum of the absolute value of the difference between atomic polarizabilities of all bonded atoms in the molecule (including implicit hydrogens) with polarizabilities, a\_acc = hydrogen bond acceptor atoms (not counting acidic atoms but counting atoms that are both hydrogen bond donors and acceptors such as –OH), a\_heavy = number of heavy atoms, logP(*o*/*w*) = Log of the octanol/water partition coefficient (including implicit hydrogens), logS = Log of the aqueous solubility (mol/L), TPSA = Polar surface area (Å<sup>2</sup> ) calculated using group contributions to approximate the polar surface area from connection table information only, Weight = molecular weight (including implicit hydrogens) in atomic mass units with atomic weights, wienerPol = wiener polarity number (half the sum of distance matrix entries with a value of 3), r<sup>2</sup> = Correlation coefficient, RMSE = Root mean square error, q<sup>2</sup> = Cross-validated r<sup>2</sup> , RMSE(LOO) = Cross-validated RMSE.

#### *2.2. Structure-Based Screening of Filtered Compounds against α-Glucosidase*

Previously we have generated the 3D-coordinates of *S. cerevisiae α*-glucosidase by homology modeling [17–19], which is composed of 579 residues. According to Ramachandran plot, the model possesses excellent stereochemical properties. Out of 579, 444 (86.7%), 63 (12.3%) and 3 (0.6%) residues lied in the most favored, additionally allowed and generously allowed regions, respectively. While two residues (0.4%) (Ala278 and Thr566) are present in disallowed regions which are not a part of active site. ERRAT showed 93.52 quality factor of the model, while in verify3D plot, 95.5% residues showed average 3D-1D score of 0.2. The stereochemical and geometric properties of the model is good and can be used in the structure-based filtration of compounds. The structural geometry and topology of the model was found to be similar to the structural topology of its template (Figure 2). The catalytic residues are conserved among *S. cerevisiae* isomaltase and *S. cerevisiae α*-glucosidase. The sequence alignment of the model and the template is shown in (Figure S1, Supporting Information). The substrate molecule (isomaltose taken from PDB ID 3AXH [22]) was manually docked to deduce the important catalytic residues of protein. Asp214, Glu276,

and Asp349 constitutes the catalytic triad for substrate catalysis, where Asp214 and Glu276 serves as nucleophile and proton donor, respectively and Asp349 stabilized the transition state of substrate molecule. The lining of the active site is composed to Asp68, Tyr71, Val108, His111, Phe157, Phe158, Phe177, Gln181, Arg212, Thr215, Ala278, Phe300, Arg312, His348, Gln350, Asp408, Arg439, and Arg443 that mediates multiple hydrophilic and hydrophobic interactions with the substrate and the competitive inhibitor of α-glucosidase 'acarbose'. Moreover, we have recognized that ten water molecules (1021,1026, 1056, 1058, 1061, 1087, 1102, 1122, 1174, and 1228) in the active site play important role in proteinsubstrate/protein-inhibitor bridging (Figure 3).

**Figure 2.** The superimposed view of model (coral ribbons) and templates [3AXH (golden ribbon) and 3A47 (Cyan ribbon)] shows the structural topology of the model is similar to its templates. The active site residues are shown in stick model. The substrate molecule (isomaltose) is depicted in green stick model.

For structure-based screening, ZINC database [22] was filtered according to the physicochemical properties of the substrate molecule (isomaltose) and 6609 compounds were matched with the given parameters. A dataset of 6609 compounds with 38 known inhibitors (αGIs used in QSAR modeling) was docked in the active site of *S. cerevisiae α*-glucosidase to determine their binding potential with α-glucosidase enzyme. The known inhibitors (αGIs) were embedded in the screening library to test the screening accuracy of docking method. The virtual screening accuracy was tested by using two metrics namely enrichment factor (EF) and Receiver operating characteristic-curve (ROC-curve). The analysis of EF and ROC-curve is discussed in Supporting Information. MOE showed >7, >47 and >78%EF in top 1%, 5% and 10% screened library. Additionally, ROC-curve shows area under the curve (AUC) value of 0.94, reflecting good virtual screening performance of MOE. The ROC-curve and the correlation of inhibitory activities of known inhibitors with their docking score is displayed in Figures S2 and S3, Supporting Information.

**Figure 3.** The binding mode of substrate is shown in the active site of enzyme. The binding interactions are highlighted in box. Hydrogen bonds are shown in green lines.

> The docked library was ranked according to the docking score and top 10% of the screened compounds (>600) were selected for post-docking filtration. Based on docking ranks, scores and calculated protein-ligand interaction fingerprints (PLIF), 202 compounds were considered 'best' and their pharmacokinetics (absorption, distribution, metabolism, excretion and toxicity) profiling was carried out by computational tools [SwissADME (http://www.swissadme.ch/ accessed on 15 December 2020) and admetSAR (http://lmmd. ecust.edu.cn/admetsar2 accessed on 15 December 2020)].

#### *2.3. ADMET Analysis of Selected Hits*

The top ranked docked compounds (202 best hits) were selected for their ADMET prediction which showed that 142/202 compounds exhibited acceptable pharmacokinetic properties (named as compounds **1**–**142**, Table S1, Supporting Information). Those compounds do not penetrate blood brain barrier (BBB), display high human intestinal absorption (HIA), revealed no AMES toxicity and carcinogenicity. The calculated oral toxicity showed that these compounds fall in category III of acute oral toxicity (LD<sup>50</sup> = >500 mg/kg to <5000 mg/kg), indicating that the compounds are non-toxic. The predicted RAT acute toxicity reflect that the compounds show low in vivo toxicity up to

the dose of ≥2 mol/kg. The synthetic accessibility (SA) values of these compounds are also ≤5 indicating that these compounds are easily synthesizable.

### *2.4. Protein-Ligand Interaction Analysis of 142 Compounds*

According to ADMET profile, compounds **1**–**142** were retrieved as good inhibitors. The modes of binding of all those (142) compounds were analyzed in the active site of α-glucosidase to select the potential inhibitors of α-glucosidase. The compounds were sorted on the basis of higher number of hydrogen bonding pattern within the active site, and compounds with number of hydrogen bonds ≥ 3, were selected. Subsequently eight molecules (**20, 28, 48, 63, 94, 112, 135 and 140**) were retrieved as 'Hits' which specifically interacted with one or two residues of catalytic triad (Asp214-Glu276-Asp349). Therefore, based on docking scores and binding interactions, those 'Hits' were considered as most active inhibitors. The predicted protein-ligand interactions of 142 compounds are shown in Figure 4. The docking ranks, scores, and molecular interactions of each compound are summarized in Table S2 (Supporting Information).


**Figure 4.** The predicted interaction pattern of compounds **1**–**142.** The binding interactions are shown in 3-color-scale (green through red colour scheme) where red and green indicates least and high number of bonds, respectively, HB = Hydrogen bonds, II = ionic interactions, WB = water bridging.

#### *2.5. The Binding Potential of High Active Hits*

The analysis of binding interactions reflect that eight compounds (**20, 28, 48, 63, 94, 112, 135** and **140**) can act as potential inhibitors of α-glucosidase enzyme. Those compounds exhibited good binding interactions within the active site of enzyme. The docked view of compound **20** showed that Glu276 and Asp214 of catalytic triad (CT) and Thr215 formed hydrogen bonds with the amino group of the compound. Moreover, the side chains of Asp68, Asp214 and Glu276 mediated ionic interaction with the amino and the nitrile groups of the compound. Additionally, a water molecule was involved in the proteinligand bridging. These interactions are responsible for the higher binding affinity of the compound [docking score (DS) = −18.35 kcal/mol and docking rank (DR) = 34). The triazolopyrimidine ring of compound **28** (DS= −18.04 kcal/mol, DR = 46) interacted with the side chain of Asp349 of CT, while the hydroxyl group of the compound formed H-bonds with the side chains of Asp214 and Arg212. Moreover, Asp349 and Phe157 provide ionic and π-π interactions to the compound, respectively. The morpholine ring of **48** (DR = 77, and DS = −17.50 kcal/mol) interacted with Asp214 and Arg439 through H-bonds, while Asp214 also mediated ionic interaction with the compound. The amino group of the compound **63** (DR = 99, and DS = −16.99 kcal/mol) formed H-bonds with the side chains of Asp214 and Thr215, while the sulfone moiety of **63** interacted with the amino group of Thr215 via H-bond. The docked view of compound **94** (DS = 16.47 kcal/mol, DR = 143) showed that the thioxoimidazolidinone moeity of **94** interacts with multiple residues including Asp68, Asp214 and His111 through H-bonds. Additionally, two water molecules also stabilize the compound by H-bonding. The binding mode of compound **112** showed that the triazolopyridine ring of the compound mediated ionic interactions with Asp349, while the substituted hydroxyl group of the compound formed H-bonds with Glu276 and Arg212. The compound depicted DS and DR of −16.27 kcal/mol and 167, respectively. The docked view of **135** showed that one of the nitrile group formed H-bond with His348, while amino group and piperidine nitrogen of the compound mediated H-bonds Asp68 and Glu276, respectively. Moreover, the side chains of Asp68 and Glu276 also offered ionic interaction to the compound. The compound exhibited DS = −16.02 kcal/mol and DR =194. The binding mode of compound **140** (DS = −15.95 kcal/mol, DR = 200) depicted that the hydroxytetrahydropyrimidinone ring of **140** binds with Asp214, Asp349, Arg212 and His348 through H-bonds.

The binding interactions reflect that these compounds possess strong binding potential with α-glucosidase. The chemical structures of the selected hits, and their binding interactions with the active site residues are shown in Table 2 and docking results are tabulated in Table S2. The 3D-structures of selected hits are given in Table S3. The binding modes of the selected hits are depicted in Figure 5. Therefore, by using virtual screening, we selected most appropriate binders of α-Glucosidase enzyme. Thus, this study will be useful in the identification of novel and more potent anti-diabetic compounds.


**Table 2.** The Chemical Structures and Binding Interactions of Eight Hits.


**Table 2.** *Cont.*

HB= Hydrogen bonds, II = Ionic interactions, WB = water bridging.

**Figure 5.** The binding modes of compounds **20, 28, 48, 63, 94, 112, 135** and **140** are depicted in the active site of α-glucosidase. Ligands are shown in green stick model, H-bonds are displayed in black lines, the active site residues are presented in orange stick model.

#### *2.6. Prediction of α-Glucosidase Inhibitory Activities of Compounds (20, 28, 48, 63, 94, 112, 135 and 140) by QSAR Model*

The ADMET profile, and predicted binding interactions of compounds 1-142, led us to identify eight new compounds as good inhibitors of α-glucosidase, therefore their α-glucosidase inhibitory activities were predicted by QSAR model. Prior to the prediction, the model was validated by a set of thirty-two compounds (regarded as Test set 2). The compounds in Test set 2 were randomly selected from literature with α-glucosidase inhibitory activities. The developed QSAR model predicted their activities with r<sup>2</sup> value of 0.62, indicating that the QSAR model can accurately predict the α-glucosidase inhibitory potential of compounds. The results are shown in Table S4 and Figure S4.

The predicted activities (pIC<sup>50</sup> values) of compounds **20, 28, 48, 63, 94, 112, 135** and **140** are 5.57, 6.79, 4.04, 4.53, 4.71, 5.85, 5.39, and 4.72, respectively. It reflects that these compounds will serve as better inhibitors when tested in vitro.

#### **3. Methods and Materials**

#### *3.1. Selection of Data Set for QSAR Analysis*

A set of thirty-eight structurally diverse compounds (α**GIs**) with α-glucosidase inhibitory activities was retrieved from literature [42–45]. α**GIs** covers a good range of inhibitory activities. α**GI1**- α**GI16**, α**GI17**- α**GI33**, and α**GI34**- α**GI38**, are reported by Wang et al. 2017, Taha et al. 2008, and Alhassan et al. 2018, respectively [42–45]. Based on the diversity of structure and the activity, the data was segregated into training and test sets of 32 and 6 compounds, respectively. The training set was used to develop QSAR model whereas test set was used for the validation of the generated model. The range of biological activities and chemical diversity of compounds was evenly spanned in both sets. The IC<sup>50</sup> values were converted into negative logarithmic scale (pIC50= −logIC50) which cover an interval of 3 log units (Table 1). The 3D-cordinates of compounds were built by MOE [46]. Partial charges were applied on compounds and their structures were minimized by MMFF94 force field until gradient was reached to 0.1 kcal/mol/Å<sup>2</sup> [46].

#### The Generation and Validation of QSAR Model

QSAR modeling was performed on MOE. Initially 192 2D-descriptors were calculated for each compound of training set. The best descriptors were selected by using MOE "QuaSAR-Contingency" applications. Descriptors with contingency coefficient above 0.6 and Cramer's, uncertainty, and correlation coefficients above 0.2 were selected for QSAR which proposed nine 2D-descriptors including apol (Sum of the atomic polarizabilities), bpol (Sum of the absolute value of the difference between atomic polarizabilities of all bonded atoms in the molecule), a\_acc (Number of hydrogen bond acceptor atoms), a\_heavy (Number of heavy atoms), logP(*o*/*w*) (Log of the octanol/water partition coefficient), logS (Log of the aqueous solubility (mol/L)), TPSA (topological polar surface area (Å2) c), Weight (molecular weight) and wienerPol (Wiener polarity number) as "best" for our dataset. These descriptors are widely used to predict biological activity and ADMET properties of dataset [46]. QSAR model was built by selecting the inhibitory activities of compounds (as descriptor variable) and calculated descriptors as model fields. Regression analysis was performed on the training set, and r<sup>2</sup> and RMSE values of the fit were obtained. The generated QSAR model was validated by cross validation by leave-one-out method. QSAR fit was used to validate the model by cross validation method. The predictive activities and the residual values of each compound of training set was assessed. Residual value, Z-score and predicted values were calculated for validation and cross-validation of model. The results were interpreted by drawing correlation plot between predicted values (actual activity) on x-axis versus predicted IC<sup>50</sup> values on y-axis. The developed model was used to predict the biological activities of test set and test set 2 compounds in terms of correlation coefficient (r<sup>2</sup> ) between the experimental and predicted activities of test set. The developed QSAR model was used for the prediction of biological activities of selected Hits.

#### *3.2. Filtration of ZINC Database for Virtual Screening*

A dataset from ZINC database [47] was selected for virtual screening against αglucosidase. Compounds were searched by using filter parameters calculated by the physicochemical properties of substrate molecule (molecular weight = 150–500, range of xlogP= −4 to 5, number of rotatable bonds = 0–8, maximum topological surface area = 150, number of hydrogen bond donor and acceptor = 0–10, polar desolvation = −400 to 1 kcal/mol, apolar desolvation = −100 to 40 kcal/mol and net charge = −4 to 5). 6609 compounds were matched with the filter parameters that were selected and imported into MOE database, their protonation states were set according to neutral pH, partial charges were applied, and energy was minimized (as described above).

#### *3.3. Docking Based Screening*

The crystallographic structure of *Saccharomyces cerevisiae* α-glucosidase is unavailable. Thus, α-glucosidase model was generated in our previous studies [17–19]. The homology model was used in the screening of 6609 compounds (collected from ZINC database) and 38 known inhibitors by MOE docking suit [46]. The protonation state of protein model was assigned according to neutral pH, partial charges were applied, and 3D-coordinates of the model was minimized by MMFF94x force field by retaining all heavy atoms fixed until RMSD gradient was reached to 0.1 kcal/mol/Å<sup>2</sup> . Triangle matcher docking algorithm was used with London dG scoring function and GBVI/WSA rescoring method. Force-field based scoring function was applied for post-docking refinement. Later, protein-ligand interaction fingerprint (PLIF) of MOE was used to calculate the binding interactions of compounds within the active site of α-glucosidase. PLIF calculates hydrogen bonding, ionic and hydrophobic interactions [46].

#### *3.4. Pharmacokinetic (ADMET) Analysis*

Pharmacokinetic properties of compounds were calculated by SwissADME (http: //www.swissadme.ch/ (accessed on 15 December 2020)) and admetSAR servers (http: //lmmd.ecust.edu.cn/ (accessed on 15 December 2020)). These servers are widely used to calculate the absorption, distribution, metabolism, excretion and toxicity (ADMET) of compounds. They predict the drug likeness, lead likeness and synthetic feasibility of compounds. These servers predict the ADMET/drug likeness/lead likeness/pharmacokinetics of compounds via different similarity searching methods or ADMET-QSAR models which are developed by the structures of several drugs or drug like compounds (curated from literature). For ADMET calculation, SMILE format of each compound was uploaded on servers and the obtained results were saved in tabular form.

#### **4. Conclusions**

α-Glucosidase inhibitors are class of oral medications which have a promising role in the management of glycemic control in T2DM. In the present study, computational techniques including QSAR modeling and structure-based virtual screening were used to identify potent inhibitors of α-glucosidase. The developed linear QSAR model displayed high predictability with correlation coefficient value (r<sup>2</sup> = 0.88). The test set validated the model with significant r<sup>2</sup> (0.82) values. Additionally, virtual screening of 6609 compounds was performed by molecular docking on *S. cerevisiae* α-glucosidase homology model. The best binders were further screened by ADMET profiling and interaction analysis. Eventually eight compounds were selected with high binding affinity for α-glucosidase. Results indicates that those in-silico identified α-glucosidase inhibitors can block the biological activity of α-glucosidase when tested in vitro.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/ph14050482/s1, Procedure of Homology Modeling and Calculation of Enrichment Factor (EF) and Receiver operating characteristic-curve (ROC-curve), Figure S1: The sequence alignment of S. cerevisiae α-glucosidase model and templates, Figure S2: Receiver Operating Curve (ROC) for α-glucosidase obtained from MOE, Figure S3: correlation of inhibitory activities of known inhibitors with their docking score, Table S1: ADMET Properties of Selected 142 Hits, Table S2: Docking scores, docking ranks and protein-ligand interactions of 142 Hits, Table S3: 3D-structures of selected hits, Table S4: The predicted activities of thirty-two compounds (Test set 2), Figure S4: The correlation plot of Test set 2.

**Author Contributions:** Conceptualization, S.A.H. and A.A.-H.; Methodology, S.A.H. and S.J.; Software, S.J.; Formal Analysis, A.K.; Writing—Original Draft Preparation, S.A.H. and A.K.; Writing— Review & Editing, S.A.H. and A.K.; Supervision, A.A.-H.; Funding Acquisition, A.A.-H. All authors have read and agreed to the published version of the manuscript.

**Funding:** The project was supported by grant from The Oman Research Council (TRC) through the funded project BFP/RGP/HSS/18/018).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data is contained within the article or Supplementary Materials.

**Acknowledgments:** The authors would like to thank the University of Nizwa for the generous support of this project.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Marine-Derived Natural Products as ATP-Competitive mTOR Kinase Inhibitors for Cancer Therapeutics**

**Shraddha Parate 1 , Vikas Kumar 2 , Gihwan Lee 1 , Shailima Rampogu 2 , Jong Chan Hong 1, \* and Keun Woo Lee 2, \***


**Abstract:** The mammalian target of rapamycin (mTOR) is a serine/threonine kinase portraying a quintessential role in cellular proliferation and survival. Aberrations in the mTOR signaling pathway have been reported in numerous cancers including thyroid, lung, gastric and ovarian cancer, thus making it a therapeutic target. To attain this objective, an in silico investigation was designed, employing a pharmacophore modeling approach. A structure-based pharmacophore (SBP) model exploiting the key features of a selective mTOR inhibitor, Torkinib directed at the ATP-binding pocket was generated. A Marine Natural Products (MNP) library was screened using SBP model as a query. The retrieved compounds after consequent drug-likeness filtration were subjected to molecular docking with mTOR, thus revealing four MNPs with better scores than Torkinib. Successive refinement via molecular dynamics simulations demonstrated that the hits formed crucial interactions with key residues of the pocket. Furthermore, the four identified hits exhibited good binding free energy scores through MM-PBSA calculations and the subsequent in silico toxicity assessments displayed three hits deemed essentially non-carcinogenic and non-mutagenic. The hits presented in this investigation could act as potent ATP-competitive mTOR inhibitors, representing a platform for the future discovery of drugs from marine natural origin.

**Keywords:** mTOR kinase; marine natural products; ATP-competitive inhibitors; structure-based pharmacophore modeling; virtual screening; molecular docking; molecular dynamics simulations; binding free energy; in silico ADMET

#### **1. Introduction**

The growth factors, nutrients and energy levels in cells are key determining factors of cellular growth and proliferation [1]. Dysregulation of the phosphoinositide 3-kinase (PI3K)/AKT/mammalian target of rapamycin (mTOR) axis has been heavily implicated in tumorigenesis and the progression of numerous cancers, including lung, thyroid, ovarian and gastric cancer [2–5]. Therefore, targeting mTOR signaling represents an attractive therapy in cancer. mTOR is a serine/threonine kinase (UniProt ID: P42345), first discovered in budding yeast Saccharomyces cerevisiae [6] and functions as a cardinal regulator of cell growth, proliferation, metabolism, energy homeostasis, angiogenesis and survival [7,8]. In mammalian cells, mTOR exists in two evolutionarily conserved complexes: mTORC1, which regulates protein synthesis through the phosphorylation of p70S6K1, and 4E-BP1 and mTORC2, which regulate cell survival and proliferation through the phosphorylation of AKT/PKB [9,10]. The mTORC1 comprises regulatory-associated protein of mTOR (RAPTOR), lethal with SEC13 protein 8 (LST8), proline-rich substrate of 40 kDa (PRAS40) and domain-containing mTOR-interacting protein (DEPTOR), while mTORC2 consists

**Citation:** Parate, S.; Kumar, V.; Lee, G.; Rampogu, S.; Hong, J.C.; Lee, K.W. Marine-Derived Natural Products as ATP-Competitive mTOR Kinase Inhibitors for Cancer Therapeutics. *Pharmaceuticals* **2021**, *14*, 282. https://doi.org/10.3390/ph14030282

Academic Editor: Osvaldo Andrade Santos-Filho

Received: 15 February 2021 Accepted: 18 March 2021 Published: 21 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of rapamycin-insensitive companion of mTOR (RICTOR), LST8, stress-activated protein kinase interacting protein 1 (SIN1), DEPTOR and protein observed with RICTOR (PRO-TOR) [6]. Deregulation of mTORC1/ mTORC2 both upstream and downstream is implicated in various cancers, including breast, ovarian, prostate and lung cancer [11].

Currently, numerous pharmacological possibilities have been developed to inhibit mTOR, resulting in three generations of mTOR inhibitors. Rapamycin and its analogs (rapalogs) belong to the first generation of mTORC1 inhibitors, approved for different cancer treatments [12]. In spite of being approved, these inhibitors cause the stabilization of the disease and not tumor regression, thereby behaving as cytostatic and not as cytotoxic. Additionally, their continued use results in copious adverse effects causing suppression of the immune system, reduction in male fertility and hematological toxicity [13]. The second generation of mTOR inhibitors act as ATP-competitive inhibitors of the catalytic kinase domain controlling the activity of both mTORC1 and mTORC2 [14]. Various resistant mutations interfering with mTOR drug binding are observed in rapalogs as well as ATP-competitive kinase domain inhibitors. The third generation of mTOR inhibitors consist of bivalent drugs binding simultaneously to the regulator and mTOR catalytic sites thereby supplementing to eliminate the resistant mutations [15,16]. In recent years, several second generation of mTOR inhibitors were identified as selective ATP-competitive kinase domain inhibitors including INK-128, OSI-027 and CC-223. Moreover, dual mTOR/PI3K inhibitors, including BEZ235, PF-04691502 and GSK2126458, were also discovered as effective inhibitors. Nevertheless, these inhibitors exhibit detrimental effects including thrombocytopenia, depression, weight loss, skin rash and mucositis [17,18]. Hence, there is an emergent need to discover potent mTOR kinase domain inhibitors as therapeutic candidates for cancerous ailments.

Natural products are predominant sources of active ingredients in medicine, demonstrating varying structural diversity and exhibiting greater potential druggable pharmacophores than synthetic molecules [19,20]. Moreover, lead compounds of natural origin have proven to show significant inhibitory activities on mTOR kinase domain in the past and this provided a structural basis for identifying natural products as potent mTOR inhibitors [21–23]. Natural products from marine resources have recently shown therapeutic potential in the development of anti-cancer drugs [24–30]. It is noteworthy to mention that 7 compounds from marine organisms have been approved as pharmaceuticals for marketing, 23 compounds are in various phases of clinical trials and about 1000 compounds are undergoing preclinical studies [25]. Four marine compounds, namely cytarabine (Cytosars), trabectedin (Yondeliss), eribulin mesylate (Halavens) and the conjugated antibody brentuximab vedotin (Acentriss), are currently utilized as anti-cancer therapeutics [25]. Marine natural extracts have also exhibited significant anti-cancer effects in recent years [24,31–34].

The aforementioned perspectives prompted us to identify natural compounds from marine resources as potential therapeutics targeted to mTOR for treatment in cancer. In accordance with it, we have carried out an in silico study to identify mTOR inhibitors via structure-based pharmacophore modeling approach. We have used pharmacophorebased virtual screening followed by molecular docking to select candidates from a Marine Natural Product (MNP) library. Subsequently, the drug-like marine compounds showing the best docking scores and molecular interactions with the kinase domain of mTOR were further refined by molecular dynamics (MD) simulations and Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA) analysis. The compounds demonstrating good binding affinity scores, as revealed by MM-PBSA were confirmed as final hits and reported in this study as potential ATP-competitive mTOR kinase domain inhibitors for cancer therapeutic treatment.

#### **2. Results**

#### *2.1. Structure-Based Pharmacophore Model*

A structure-based pharmacophore was generated from the crystallographic structure of mTOR kinase complexed with PP242 inhibitor, in which the key features of PP242 binding

with mTOR were exploited. Accordingly, a total of six pharmacophore models were generated with hydrogen bond donor (HBD), hydrogen bond acceptor (HBA) and hydrophobic (Hy) as common indispensable features (Table 1). The *Pharmacophore\_01* with five features (1 HBA, 2 HBD and 2 Hy) and with the highest selectivity score of 9.2973 was observed to be the best pharmacophore model among the six generated pharmacophores and hence selected for further analysis. Upon meticulous examination, *Pharmacophore\_01* displayed requisite features complementing the key residues for binding—Asp2195, Gly2238, Val2240 and Ile2356. The HBA feature complements with the essential residue Val2240, where PP242 is responsible for forming a hydrogen bond. The pyrazolopyrimidine scaffold of PP242 also maps on to the two HBD features by hydrogen bonding with two vital residues— Asp2195 and Gly2238. A previously published study suggests HBA and HBD as crucial features required for mTOR inhibition [21]. Moreover, the two Hy features complement the PP242 binding with Ile2356 residue via π-alkyl interactions. The Hy pharmacophoric feature was also reported as an essential feature in an earlier study on nanomolar mTOR inhibitors [35]. Therefore, *Pharmacophore\_01* was escalated for further analysis (Figure 1). Pharmaceuticals 2021, 14, x FOR PEER REVIEW 3 of 19 2. Results 2.1. Structure-Based Pharmacophore Model A structure-based pharmacophore was generated from the crystallographic structure of mTOR kinase complexed with PP242 inhibitor, in which the key features of PP242 binding with mTOR were exploited. Accordingly, a total of six pharmacophore models were generated with hydrogen bond donor (HBD), hydrogen bond acceptor (HBA) and hydrophobic (Hy) as common indispensable features (Table 1). The Pharmacophore\_01 with five features (1 HBA, 2 HBD and 2 Hy) and with the highest selectivity score of 9.2973 was observed to be the best pharmacophore model among the six generated pharmacophores and hence selected for further analysis. Upon meticulous examination, Pharmacophore\_01 displayed requisite features complementing the key residues for binding—Asp2195, Gly2238, Val2240 and Ile2356. The HBA feature complements with the essential residue Val2240, where PP242 is responsible for forming a hydrogen bond. The pyrazolopyrimidine scaffold of PP242 also maps on to the two HBD features by hydrogen bonding with two vital residues—Asp2195 and Gly2238. A previously published study suggests HBA and HBD as crucial features required for mTOR inhibition [21]. Moreover, the two Hy features complement the PP242 binding with Ile2356 residue via π-alkyl interactions. The

> **Table 1.** Structure-based pharmacophore models with their generated features. on nanomolar mTOR inhibitors [35]. Therefore, Pharmacophore\_01 was escalated for fur-


Hy pharmacophoric feature was also reported as an essential feature in an earlier study

\* A: hydrogen bond acceptor (HBA); D: hydrogen bond donor (HBD); H: hydrophobic (Hy). Pharmacophore\_06 4 ADHH 6.8689 \*A: hydrogen bond acceptor (HBA); D: hydrogen bond donor (HBD); H: hydrophobic (Hy).

Figure 1. Structure-based pharmacophore model- Pharmacophore\_01. (A) Pharmacophore model generated at the catalytic site of mTOR with co-crystallized ligand, PP242. (B) Pharmacophore features mapped with the key residues of the binding pocket. (C) Interfeature distance between the mapped pharmacophore features. HBA (hydrogen bond acceptor); HBD (hydrogen bond donor) and Hy (hydrophobic). **Figure 1.** Structure-based pharmacophore model- *Pharmacophore\_01*. (**A**) Pharmacophore model generated at the catalytic site of mTOR with co-crystallized ligand, PP242. (**B**) Pharmacophore features mapped with the key residues of the binding pocket. (**C**) Interfeature distance between the mapped pharmacophore features. HBA (hydrogen bond acceptor); HBD (hydrogen bond donor) and Hy (hydrophobic).

#### *2.2. Decoy Set Validation of the Structure-Based Pharmacophore Model*

The selected model *Pharmacophore\_01* was evaluated for its efficiency in retrieving active mTOR compounds from a given database of active and inactive molecules. This validation was prompted by screening an external database (D) of 300 compounds, with 50 active compounds (A). With 61 hits retrieved from the database (Ht), active compounds obtained were 49 (Ha). The goodness of fit (GF) score was calculated as 0.80, thereby confirming that *Pharmacophore\_01* can predict active compounds from a given dataset reasonably well (Table 2).


**Table 2.** Decoy set validation of *Pharmacophore\_01* from an external database composed of active and inactive mTOR inhibitors.

#### *2.3. Drug-Like Marine Compounds Retrieved by Virtual Screening*

The validated model *Pharmacophore\_01* mapped 3019 compounds from the Marine Natural Products (MNP) library of 14,492 compounds. Subsequent filtering of mapped compounds by Lipinski's Rule of Five (Ro5) and absorption, distribution, metabolism, excretion and toxicity (ADMET) properties led to further reducing the amount to 135 compounds. These 135 marine drug-like compounds along with reference inhibitor PP242 were escalated for molecular docking with mTOR crystallographic structure and their interactions with residues Leu2185, Asp2195, Ile2237, Gly2238, Trp2239, Val2240, Thr2245, Met2345, Pharmaceuticals Leu2354, Ile2356, and Asp2357 were scrutinized (Figure 2). <sup>2021</sup>, 14, x FOR PEER REVIEW 5 of 19

Figure 2. Illustration of the stages involved in the retrieval of potential drug-like compounds from Marine Natural Products (MNP) library using the structure-based pharmacophore model. 2.4. Molecular Docking of Retrieved Marine Drug-Like Compounds with mTOR Kinase **Figure 2.** Illustration of the stages involved in the retrieval of potential drug-like compounds from Marine Natural Products (MNP) library using the structure-based pharmacophore model.

itor. The reference compound displayed a lower Gold score of 63.20 as compared with marine compounds hereafter referred to as MNP1, MNP2, MNP3 and MNP4 exhibiting

The molecular docking process demarcates on the binding affinity, mode of ligand binding in the target protein pocket and also elucidates the interactions of compounds

#### *2.4. Molecular Docking of Retrieved Marine Drug-Like Compounds with mTOR Kinase*

The molecular docking process demarcates on the binding affinity, mode of ligand binding in the target protein pocket and also elucidates the interactions of compounds with essential residues. The performance of GOLD software was evaluated by re-docking the co-crystallized ligand PP242 into mTOR binding pocket, resulting in an acceptable RMSD of 0.71 Å (Figure S1). Docking of screened 135 marine drug-like compounds with mTOR kinase domain was then carried out along with reference inhibitor PP242. A total of four marine compounds demonstrated higher Gold scores than PP242 reference inhibitor. The reference compound displayed a lower Gold score of 63.20 as compared with marine compounds hereafter referred to as MNP1, MNP2, MNP3 and MNP4 exhibiting higher Gold scores of 65.48, 65.41, 64.72 and 63.75, respectively (Table 3). The four compounds also demonstrated interactions with the aforementioned residues of the ATP-binding pocket of mTOR kinase domain. A total of 15 molecules exhibited lower docking scores than Torkinib and, therefore, were not considered for further evaluation. However, it was observed that these 15 marine molecules also target similar residues of the mTOR pocket, as seen for the above four compounds (Table S1). Therefore, these four marine compounds and the reference PP242 inhibitor were taken forward for MD simulations to confirm on their stabilities with the mTOR binding pocket.

**Table 3.** The docking scores and intermolecular interactions of reference PP242 and Marine Natural Product (MNP) library compounds with mTOR kinase domain (PDB ID: 4JT5).


\* CAS: Chemical Abstracts Service.

#### *2.5. Binding Mode and Binding Free Energy Analysis of Identified Marine Compounds by Molecular Dynamics Simulations*

MD simulations employing GROningen Machine for Chemical Simulations (GRO-MACS) package were executed for the four identified marine compounds docked with the mTOR kinase to comprehend the dynamic behavior at the atomistic level. Simulations were also supplemented with the calculation of binding free energies for the identified hits by MM-PBSA methodology and the compounds were ranked as hits accordingly (Table S2). To gain insight into the binding mode of identified compounds inside the mTOR ATP-binding pocket, the representative structures from the last 5 ns of stable MD trajectories were extracted and subsequently superimposed. Upon scrupulous analysis, it was observed that the identified compounds exhibited a similar binding mode as the PP242 co-crystallized ligand of mTOR (Figure 3).

Pharmaceuticals 2021, 14, x FOR PEER REVIEW 7 of 19

Figure 3. Binding mode of reference PP242 and identified Marine Natural Product (MNP) library hits within the ATP-binding pocket of mTOR kinase domain. **Figure 3.** Binding mode of reference PP242 and identified Marine Natural Product (MNP) library hits within the ATP-binding pocket of mTOR kinase domain.

#### 2.6. Characteristic Binding Interaction and Binding Free Energy Analysis of the Confirmed Marine Hits with mTOR ATP-Binding Pocket Residues *2.6. Characteristic Binding Interaction and Binding Free Energy Analysis of the Confirmed Marine Hits with mTOR ATP-Binding Pocket Residues*

#### 2.6.1. mTOR-Hit1 Interaction 2.6.1. mTOR-Hit1 Interaction

The marine compound MNP2 acquired from the docking analysis (Table 3) exhibited the highest BFE of −101.187+/−17.842 kJ/mol, as investigated by MM-PBSA calculations and, therefore, referred to as Hit1 (Table S2). The estimated BFE gives insight on the diverse components of interaction energy contributing to Hit1 binding. Both electrostatic and van der Waals components contribute a major role in the binding of Hit1 with the ATP-binding pocket of mTOR, where the van der Waals contribution (−151.961+/−11.779 kJ/mol) is higher than the electrostatic component (−101.485+/−16.574 kJ/mol). The solvent accessible surface area (SASA) provides a slightly favorable contribution towards the binding of Hit1 with mTOR (−17.534+/−0.909 kJ/mol). Energy decomposition analysis led to the identification of vital residues contributing to the binding of Hit1 with mTOR (Figure 4A and 5A, Table S2). It was observed that the major contribution for Hit1 binding was from van der Waals interaction with residue Trp2239 (−8.5 kJ/mol) and hydrophobic interaction with residue Ile2356 (−7.9 kJ/mol) which was consistent with its binding mode (Figure 5B). The binding of Hit1 in the mTOR ATP-binding pocket is rendered by hydrogen bond interactions (Figures 4A and 5B) with Asp2195 (bond length: 1.72 Å) and Asp2357 (bond length: 1.96 Å). Additionally, the benzene ring of Hit1 interacted with key residues Tyr2225 (bond length: 4.25 Å), Ile2237 (bond length: 4.70 Å), Val2240 (bond length: 4.94 Å) and Ile2356 (bond length: 4.47 Å) via π-alkyl hydrophobic interactions. Moreover, mTOR residues Ile2163, Leu2185, Lys2187, Val2198, Gly2238, Trp2239, Met2345 and Phe2358 form van der Waals interactions with Hit1 (Figure 5B). The marine compound MNP2 acquired from the docking analysis (Table 3) exhibited the highest BFE of −101.187 ± 17.842 kJ/mol, as investigated by MM-PBSA calculations and, therefore, referred to as Hit1 (Table S2). The estimated BFE gives insight on the diverse components of interaction energy contributing to Hit1 binding. Both electrostatic and van der Waals components contribute a major role in the binding of Hit1 with the ATP-binding pocket of mTOR, where the van der Waals contribution (−151.961 ± 11.779 kJ/mol) is higher than the electrostatic component (−101.485 ± 16.574 kJ/mol). The solvent accessible surface area (SASA) provides a slightly favorable contribution towards the binding of Hit1 with mTOR (−17.534 ± 0.909 kJ/mol). Energy decomposition analysis led to the identification of vital residues contributing to the binding of Hit1 with mTOR (Figures 4A and 5A, Table S2). It was observed that the major contribution for Hit1 binding was from van der Waals interaction with residue Trp2239 (−8.5 kJ/mol) and hydrophobic interaction with residue Ile2356 (−7.9 kJ/mol) which was consistent with its binding mode (Figure 5B). The binding of Hit1 in the mTOR ATP-binding pocket is rendered by hydrogen bond interactions (Figures 4A and 5B) with Asp2195 (bond length: 1.72 Å) and Asp2357 (bond length: 1.96 Å). Additionally, the benzene ring of Hit1 interacted with key residues Tyr2225 (bond length: 4.25 Å), Ile2237 (bond length: 4.70 Å), Val2240 (bond length: 4.94 Å) and Ile2356 (bond length: 4.47 Å) via π-alkyl hydrophobic interactions. Moreover, mTOR residues Ile2163, Leu2185, Lys2187, Val2198, Gly2238, Trp2239, Met2345 and Phe2358 form van der Waals interactions with Hit1 (Figure 5B).

#### 2.6.2. mTOR-Hit2 Interaction

The marine compound MNP3 (hereafter referred to as Hit2) attained from docking analysis (Table 3), demonstrated with BFE of −101.041 ± 20.457 kJ/mol, was observed to be

comparable with the BFE of Hit1 (Table S2). Similar to Hit1, contribution of van der Waals (−165.824 ± 15.257 kJ/mol) and electrostatic (−70.466 ± 12.892 kJ/mol) components for the binding of Hit2 with mTOR played a major role than the corresponding SASA energy. The contribution of SASA energy component (−21.107 ± 0.998 kJ/mol) for Hit2 binding was observed to be higher than that of Hit1 with mTOR ATP-binding pocket. The major residues contributing to Hit2 binding with mTOR were also observed to be similar to Hit1 contributing residues (Figure 5C), with Trp2239 (−9.3 kJ/mol) and Ile2356 (−7.0 kJ/mol) forming hydrophobic bonds with Hit2 (Figure 5D). Elucidating on the binding interaction of Hit2 predicted by MD analysis, it was noticed that Hit2 formed three hydrogen bonds (Figures 4B and 5D) with residues Lys2187 (bond length: 1.75 Å), Thr2245 (bond length: 1.80 Å) and Asp2357 (bond length: 2.45 Å). Moreover, the bromo-phenoxy group interacts with residues Leu2185 (bond length: 4.50 Å), Ile2237 (bond length: 3.77 Å) and Trp2239 (bond length: 5.45 Å) via alkyl and π-alkyl hydrophobic interactions. Additional interactions with residues Ala2248 (bond length: 5.18 Å) and Ile2356 (bond length: 4.37 Å) also hold Hit2 in the mTOR hydrophobic pocket via π-alkyl bonds. Furthermore, the residues Ile2163, Thr2164, Glu2190, Asp2195, Gly2238, Val2240, Cys2243, Asp2244, Met2345 and Phe2358 assist in holding Hit2 in the mTOR binding pocket firmly via van der Waals interactions (Figure 5D). Pharmaceuticals 2021, 14, x FOR PEER REVIEW 9 of 19 Waals interactions, respectively, as also observed in the above hits (Figure 5G,H). The binding interaction of Hit4 with mTOR revealed that Hit4 interacts with Asp2195 (bond length: 2.03 Å), Val2240 (bond length: 2.40 Å) and Asp2357 (bond length: 2.83 Å) via hydrogen bonds (Figures 4D and 5H). Similar interactions with Asp2195 and Val2240 via hydrogen bonds were also observed in reference structure bonding with mTOR kinase (Figure S2). Elucidating on the hydrophobic interactions of Hit4 with mTOR kinase domain residues, it was observed that Hit4 interacts with hydrophobic chamber residue Ile2163 (bond length: 4.69 Å) and inner hydrophobic pocket residue Leu2185 (bond length: 4.52 Å) via alkyl bonds. Moreover, three π-alkyl bonds (bond lengths: 4.82 Å, 5.22 Å and 5.24 Å) and one π-sigma bond (bond length: 3.73 Å) with hydrophobic chamber residue Trp2239 also contribute to Hit4-mTOR interaction. Furthermore, residues Met2199, Tyr2225, Val2227, Ile2237, Pro2241, His2242, Cys2243, Met2345, Ile2356 and Phe2358 also participate in positioning Hit4 firmly in mTOR binding pocket via van der Waals interactions (Figure 5H).

Figure 4. Interactions of (A) Hit1, (B) Hit2, (C) Hit3 and (D) Hit4 from Marine Natural Product (MNP) library with key residues of mTOR ATP-binding pocket via hydrogen bonds. The compounds and interacting residues are represented as sticks while the hydrogen bonding interactions are shown as green dashed lines. **Figure 4.** Interactions of (**A**) Hit1, (**B**) Hit2, (**C**) Hit3 and (**D**) Hit4 from Marine Natural Product (MNP) library with key residues of mTOR ATP-binding pocket via hydrogen bonds. The compounds and interacting residues are represented as sticks while the hydrogen bonding interactions are shown as green dashed lines.

#### 2.6.3. mTOR-Hit3 Interaction

The marine compound MNP1, exhibiting the highest dock score of 65.48 (Table 3) presented with the BFE of −91.924 ± 12.264 kJ/mol (Table S2) and, therefore, ranked and referred as Hit3. It was observed that the van der Waals component was responsible for majorly contributing (−171.314 ± 11.172 kJ/mol) to Hit3 binding with mTOR than electrostatic (−17.958 ± 10.086 kJ/mol) and SASA (−19.290 ± 1.055 kJ/mol) energy components. The van der Waals contribution was observed to be the highest for Hit3-mTOR interaction than mTOR interactions with other hits (Table S2). Additionally, the contribution of this component is 10-fold higher than electrostatic component (−17.958 ± 10.086 kJ/mol) for Hit3-mTOR binding. Emphasizing the major residues contributing to the binding energy of Hit3 with mTOR, van der Waals interaction via Ile2356 (-9.2 kJ/mol) and hydrophobic

interactions via Trp2239 (−11.1 kJ/mol) and Met2345 (8.9 kJ/mol) were seen to have a greater impact on binding (Figure 5E,F). Hydrogen bond analysis of the predicted binding mode demonstrated with Hit3 forming two bonds (Figures 4C and 5F) with residues Val2240 (bond length: 2.05 Å) and Asp2357 (bond length: 2.37 Å). Moreover, Hit3 interacts with Leu2185 (bond length: 5.0 Å) from the inner hydrophobic pocket and with Met2345 (bond length: 5.13 Å) of mTOR hydrophobic chamber via alkyl interactions. Two π-alkyl bonds are formed with residue Trp2239 (bond length: 4.32 Å and 4.90 Å) of the hydrophobic chamber. Besides above interactions, van der Waals bonds with residues Ile2163, Met2199, Tyr2225, Val2227, Ile2237, Gly2238, Thr2245, Leu2354, His2355, Ile2356 and Phe2358 contribute to the BFE obtained for Hit3 (Figure 5F). Pharmaceuticals 2021, 14, x FOR PEER REVIEW 10 of 19

Figure 5. Energy decomposition of individual residues from MM-PBSA contributing to total binding free energy (BFE) for (A) Hit1, (C) Hit2, (E) Hit3 and (G) Hit4, and 2D interaction of (B) Hit1, (D) Hit2, (F) Hit3 and (H) Hit4 with residues of mTOR ATP-binding pocket. The hydrogen bonding interactions are shown as green dashed lines, the hydrophobic interactions are shown as pink and purple spheres and the van der Waals interactions are displayed as light green spheres. **Figure 5.** Energy decomposition of individual residues from MM-PBSA contributing to total binding free energy (BFE) for (**A**) Hit1, (**C**) Hit2, (**E**) Hit3 and (**G**) Hit4, and 2D interaction of (**B**) Hit1, (**D**) Hit2, (**F**) Hit3 and (H) Hit4 with residues of mTOR ATP-binding pocket. The hydrogen bonding interactions are shown as green dashed lines, the hydrophobic interactions are shown as pink and purple spheres and the van der Waals interactions are displayed as light green spheres.

2.7. Evaluation of Drug-Likeness, ADME and Toxicity Properties of Identified mTOR Hits

Knowledge of the drug-likeness and ADME properties of final hits is vital, prior to

Furthermore, the Toxicity Prediction (TOPKAT) module in DS was used to evaluate the toxicity properties of identified hits. Given the structural information of a compound as a query, TOPKAT relies on the concept of Quantitative Structure–Toxicity Relationship

#### 2.6.4. mTOR-Hit4 Interaction

The marine compound MNP4 (hereafter referred to as Hit4) with lowest dock score among obtained hits (Table 3) demonstrated a BFE of −86.049 ± 14.961 kJ/mol (Table S2). Similar to other hit compounds, the van der Waals component of BFE also contributed greatly to Hit4-mTOR binding (−169.882 ± 11.839 kJ/mol). Correspondingly, the contribution of SASA energy component for Hit4 binding (−19.569 ± 0.780 kJ/mol) was in comparable range with SASA energy for Hit3 binding (−19.290 ± 1.055 kJ/mol). The interaction energy of individual residues towards Hit4-mTOR binding revealed Trp2239 (−8.5 kJ/mol) and Ile2356 (−8.1 kJ/mol) as the major contributing residues via hydrophobic and van der Waals interactions, respectively, as also observed in the above hits (Figure 5G,H). The binding interaction of Hit4 with mTOR revealed that Hit4 interacts with Asp2195 (bond length: 2.03 Å), Val2240 (bond length: 2.40 Å) and Asp2357 (bond length: 2.83 Å) via hydrogen bonds (Figures 4D and 5H). Similar interactions with Asp2195 and Val2240 via hydrogen bonds were also observed in reference structure bonding with mTOR kinase (Figure S2). Elucidating on the hydrophobic interactions of Hit4 with mTOR kinase domain residues, it was observed that Hit4 interacts with hydrophobic chamber residue Ile2163 (bond length: 4.69 Å) and inner hydrophobic pocket residue Leu2185 (bond length: 4.52 Å) via alkyl bonds. Moreover, three π-alkyl bonds (bond lengths: 4.82 Å, 5.22 Å and 5.24 Å) and one π-sigma bond (bond length: 3.73 Å) with hydrophobic chamber residue Trp2239 also contribute to Hit4-mTOR interaction. Furthermore, residues Met2199, Tyr2225, Val2227, Ile2237, Pro2241, His2242, Cys2243, Met2345, Ile2356 and Phe2358 also participate in positioning Hit4 firmly in mTOR binding pocket via van der Waals interactions (Figure 5H).

#### *2.7. Evaluation of Drug-Likeness, ADME and Toxicity Properties of Identified mTOR Hits*

Knowledge of the drug-likeness and ADME properties of final hits is vital, prior to their consideration as a lead candidate in drug development and/or their usage as anti-cancer drugs. These properties for the final hits were assessed using DS and tabulated. Furthermore, the *Toxicity Prediction (TOPKAT)* module in DS was used to evaluate the toxicity properties of identified hits. Given the structural information of a compound as a query, TOPKAT relies on the concept of Quantitative Structure–Toxicity Relationship (QSTR) models for computing toxicity properties, including AMES mutagenicity and rodent carcinogenicity, based on the National Toxicology Program (NTP) dataset. According to the U.S. NTP protocol, a compound's carcinogenicity is assessed by testing it in female and male sexes of mouse as well as rat. As per the calculations, all hits were observed to obey the Lipinski's Ro5 and displayed molecular weight of less than 500 Da with predicted octanol/water partition coefficient logP of less than 5.0 and estimated hydrogen bond donors and acceptors of less than 5 and 10, respectively (Table S2). The hits also followed the Veber's rule except for Hit2 with 10 rotatable bonds. Additionally, the identified hits satisfied the ADME properties demonstrating low blood–brain barrier (BBB) penetration, no inhibition of CYP2D6, good human intestinal absorption (HIA) as well as good aqueous solubility (Table S3). The results from TOPKAT analysis suggested that our identified hits demonstrated non-carcinogenicity in both sexes of rat models, while Hit2 appeared to be carcinogenic in mouse male (Table S4). According to the TOPKAT AMES mutagenicity analysis, Hit2 was displayed as being a mutagen while other hits were observed to be non-mutagenic (Table S4).

#### *2.8. Novelty and Source Documentation of Identified mTOR Inhibitors*

As a final assessment, the source of the identified marine hits was evaluated. For this purpose, the PubChem chemistry database (https://pubchem.ncbi.nlm.nih.gov/, (accessed on 11 February 2021)) was searched with Chemical Abstracts Service (CAS) numbers/MNP IDs of the respective hits as queries to identify if our hits have already been reported in the literature for mTOR kinase or other target proteins. From PubChem literature analysis, it was found that the compound Hit1 (7-Hydroxyceratinamine) is a cyanoformamide-containing metabolite originating from a Micronesian sponge, *Aplysinella*

sp [36]. In addition, Hit2 was found to be a marine alkaloid isolated from the sponge *Psammaplysilla purpurea* [37]. The source for Hit3 could not be identified while Hit4 was recognized as a phytolide of the *Colletotrichum boninense* fungal origin [38]. As per the literature analysis, the four identified hits have not been reported as inhibitors of mTOR kinase. Therefore, the identified hits in the present study provide valuable alternatives as therapeutic candidates for further lead optimization.

#### **3. Discussion**

Kinase proteins act as chief regulatory entities in cellular biology. Moreover, their hyperactivation leads to several pathologies, including cancer. Therefore, kinases have become essential pharmacological targets and the discovery of small molecule drugs is a predominant scientific activity to mitigate the cancerous ailments. The mammalian target of rapamycin, mTOR is a dual specificity protein kinase which phosphorylates serine/threonine as well as tyrosine residues. The deregulation of mTOR is associated with diabetes, obesity, aging and various types of cancer [39–43]. The PI3K/AKT/mTOR pathway has been largely implicated in the tumorigenesis and progression of aforementioned cancers. Specific mTOR inhibitors are currently in various stages of clinical trials [44]. The mTOR inhibitors appear to be well tolerated combined with adverse effects, including myelosuppression, metabolic abnormalities, stomatitis and skin reactions, among other abnormalities [45]. These inhibitors can be classified into first generation allosteric inhibitors such as rapamycin and its analogues [46] and second generation ATP-competitive inhibitors. Several ATP-competitive mTOR inhibitors have been discovered and being tested in clinical trials including selective inhibitors like CC-223, INK-128 and OSI-027 [47]. Unlike the former category of inhibitors, the ATP-competitive mTOR inhibitors block the ATP-binding catalytic site as well as reduce the activity of mTORC1 and mTORC2 complexes and hence has become an effective strategy to suppress both the ATP and allosteric sites. Encouraged from these efforts, we pursued our research strategy to identify natural product compounds as mTOR ATP-competitive inhibitors by applying structure-based pharmacophore modelling. Such a pharmacophore strategy works towards exploiting the key interactions between the protein residues and the bound co-crystallized ligand [48].

Marine extracts have displayed a great potential as an essential source of new drugs. Marine environments remain extensively unexplored despite being a huge source of bioactive compounds against cancer. Aquatic habitats have produced a variety of marine-derived alkaloids, triterpenoids and peptides. Intriguingly, a purple sponge extract of Turkish marine origin was shown to display promising activity against a panel of tyrosine kinases and cell lines including A549, A375, KMS-12PE and K562 cancer cell lines [32]. Similarly, the brown algae-derived polysaccharide, Fucoidan, was shown to exert anti-cancer effects not only through cell cycle arrest but also by indirectly killing cancerous cells by activation of natural killer cells and macrophages. Fucoidan was also shown to demonstrate inhibitory activity against cancer A549, MCF-7, PC-3 and SMMC-7721 cells [31]. In addition, Fascaplysin which is a bis-indole of a marine sponge demonstrated anticancer activity as CDK4 inhibitor in lung cancer cell line [46]. The antitumor potential of marine algae-derived compounds has also been extensively reviewed recently [24]. We therefore designed our study to target the mTOR dysregulation in cancer using marine-derived bioactive natural products employing a series of computational methods. Using a structurebased pharmacophore approach, we have developed a pharmacophore model from mTOR protein structure co-crystallized with ligand PP242 (Figure 1). Torkinib/PP242 is a selective ATP-competitive inhibitor of mTOR with promising anti-cancer activity over numerous cancer types [49]. A total of 14,492 compounds of marine origin were screened with the pharmacophore model as a query, deriving 3019 compounds as candidates mapping the pharmacophore model. Subsequently, a drug-like database was prepared employing Lipinski's Ro5 and ADMET rules, thereby retrieving 135 compounds (Figure 2). Molecular docking-based interaction screening of these 135 compounds resulted in the identification of four compounds with higher docking scores than PP242 and similar interactions with

the ATP-binding pocket of mTOR (Table 3). Escalating the identified four compounds to molecular dynamics simulations for observing their behavior at the atomistic level gave insights into the critical residues required for the specific mTOR inhibition.

The structure-based pharmacophore and MD analysis of the mTOR-PP242 crystal structure revealed that the inhibitor targets key residues Asp2195, Gly2238 and Val2240 via hydrogen bonds entailing essential pharmacophoric features including HBD and HBA (Figure 1B and Figure S2A). It has been elucidated in previous studies that hydrogen bonds with Asp2195 and Val2240 are indispensable for mTOR inhibitory activity [50–56]. In the current study, Hit1 was observed to retain the hydrogen bonding interaction with Asp2195 (Table 3, Figure 4A), while Hit4 interaction with both residues was preserved even after 30 ns of production simulation run (Table 3, Figure 4D). Additionally, our hits formed hydrogen bonds with catalytic hydrophilic residue Asp2357 of the mTOR ATP-binding site which offers a level of specificity for our hits towards mTOR than PI3K (Figure 4) [51–55]. The compound Hit2 also formed additional hydrogen bonds with residues Lys2187 and Thr2245 (Figure 4B) and interactions with these residues were also reported in previously published studies [22,51,52,55]. A recent study reported natural products as mTOR ATPbinding site and rapamycin binding site inhibitors derived from three databases—Marine Natural Products Library, SuperNatural II and ZINC natural products—and also provided experimental evidence against mTOR for eleven compounds [23]. However, the marine compounds identified from their studies as mTOR inhibitors are distinct from our proposed marine natural product hits. In addition to the aforementioned hydrogen bonds, our hits are also characterized by several hydrophobic and van der Waals interactions (Figure 5). In particular, interaction with residue Trp2239 of the hydrophobic chamber was observed via π-stacking bonds or van der Waals interaction for our hits, similar to Torkinib interaction with mTOR (Figure 5 and Figure S2B). Earlier studies reported that Trp2239 is not present in canonical protein kinases such as PI3K and interaction with this residue provides for mTOR inhibitor specificity over PI3K [56,57]. All of these results suggest that our hits may be selective mTOR inhibitors. From our per-residue contribution analysis, Trp2239 was observed to contribute the highest for Hit3 binding (−11.1 kJ/mol) (Figure 5E) followed by Hit2 (−9.3 kJ/mol) (Figure 5C), Hit1 (−8.53 kJ/mol) (Figure 5A) and Hit4 (−8.52 kJ/mol) (Figure 5G) binding with mTOR. The Trp2239 contributed the highest energy for binding of Hit3 among the four identified hits forming strong interactions via π hydrophobic bonds as observed (Figure 5F), thus suggesting a higher level of selectivity of Hit3 binding with mTOR. In a recently published review, the significant residues involved in ligand selectivity towards mTOR over PI3K were reviewed in detail and reported as Arg770, Glu2190 and Cys2243 [50]. Our identified Hit2 displayed bonds with Glu2190 and Cys2243 via van der Waals interactions (Figure 5D), while Cys2243 was also observed to support Hit4 via van der Waals interactions in the ATP-binding pocket of mTOR (Figure 5H). The above-mentioned analyses indicate that our identified marine hits have better binding affinity, as computed from MM-PBSA binding free energy scores, and also seem to confer comparable selectivity towards mTOR as the previously identified selective inhibitor, Torkinib. Correspondingly, the identified hits also portray the essential pharmacophoric features required for mTOR inhibition, similar to Torkinib (Figures 1 and 6).

As a final analysis, the four identified marine hits were scrutinized for their physicochemical, pharmacokinetic and toxicity properties to evaluate their in vivo disposition prior to their consideration as therapeutic anti-cancer drugs. The drug-likeness properties and ADME results demonstrated that the identified hits satisfied the Ro5 criteria, could be absorbed easily in the human intestine, show good aqueous solubility, do not inhibit the CYP2D6 enzyme and do not penetrate the BBB (Table S3). Despite presenting with good Ro5 properties, Hit2 displayed with the upper limit value (i.e., 10) in the case of Veber's rule of rotatable bonds. In addition, the hits were identified as non-carcinogenic in rodent NTP models and non-AMES mutagenic, except for Hit2, which showed carcinogenicity in the mouse male NTP model and also showed AMES mutagenicity (Table S4). Although Hit2 exhibited with a BFE of −101.041 kJ/mol (Table S2) and presented with

key intermolecular interactions with mTOR (Figures 4B and 5D), it cannot be considered as a therapeutic candidate for further drug optimization, owing to its toxicity properties. Overall, we believe the potentiality of Hit1, Hit3 and Hit4 as alternatives and their scaffolds can be further explored for developing efficient ATP-binding site mTOR inhibitors for cancer therapeutics. Though the in vitro studies of our identified marine hits are further required to clinically substantiate these findings, structure-based pharmacophore modeling and virtual screening strategy can be very useful to design potent molecules as mTOR inhibitors in future drug discovery studies. In addition, our study characterizes a platform for future discovery of novel natural chemotherapeutic drugs from marine natural habitat. Pharmaceuticals 2021, 14, x FOR PEER REVIEW 14 of 19

Figure 6. Alignment of the (A) Hit1, (B) Hit2, (C) Hit3 and (D) Hit4 with the pharmacophoric features. All hits represent the HBA (hydrogen bond acceptor), HBD (hydrogen bond donor) and Hy (hydrophobic) features of Pharmacophore\_01. **Figure 6.** Alignment of the (**A**) Hit1, (**B**) Hit2, (**C**) Hit3 and (**D**) Hit4 with the pharmacophoric features. All hits represent the HBA (hydrogen bond acceptor), HBD (hydrogen bond donor) and Hy (hydrophobic) features of *Pharmacophore\_01*.

#### 4. Materials and Methods **4. Materials and Methods**

#### 4.1. Structure-Based Pharmacophore Model Generation Receptor-based pharmacophore model delves into the catalytic site of the target pro-*4.1. Structure-Based Pharmacophore Model Generation*

tein bound with its inhibitor to identify essential pharmacophoric features effective for inhibition [58]. The structure of mTOR target protein bound with its ATP-site inhibitor PP242 (PDB ID: 4JT5, 3.45 Å) was retrieved from Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank (PDB) and considered for the generation of a structure-based pharmacophore model [56]. PP242 (also known as Torkinib) is a potent and specific inhibitor of mTOR kinase domain proven to be more effective than traditional mTOR inhibitor rapamycin [14,59]. Subsequently, the residues within 9 Å around PP242 were considered and Receptor-Ligand Pharmacophore Generation module embedded in Discovery Studio (DS) v.2018 was employed for model generation. The model with the highest selectivity score was chosen for subsequent validation. 4.2. Validation of Generated Pharmacophore Model Receptor-based pharmacophore model delves into the catalytic site of the target protein bound with its inhibitor to identify essential pharmacophoric features effective for inhibition [58]. The structure of mTOR target protein bound with its ATP-site inhibitor PP242 (PDB ID: 4JT5, 3.45 Å) was retrieved from Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank (PDB) and considered for the generation of a structure-based pharmacophore model [56]. PP242 (also known as Torkinib) is a potent and specific inhibitor of mTOR kinase domain proven to be more effective than traditional mTOR inhibitor rapamycin [14,59]. Subsequently, the residues within 9 Å around PP242 were considered and *Receptor-Ligand Pharmacophore Generation* module embedded in Discovery Studio (DS) v.2018 was employed for model generation. The model with the highest selectivity score was chosen for subsequent validation.

#### ensuring efficient retrieval of active target protein compounds. Accordingly, the model *4.2. Validation of Generated Pharmacophore Model*

chosen from the above-mentioned criteria was validated by Güner–Henry approach (decoy set) approach [60] for evaluating the robustness of the pharmacophore model on the basis of goodness of fit (GF) score in the range of 0 (null model) and 1 (ideal model) [58,61,62]. = ൬ 4<sup>൰</sup> (3 + ) × ൜1 − − − <sup>ൠ</sup> Given a particular dataset, pharmacophore validation is a quintessential criterion for ensuring efficient retrieval of active target protein compounds. Accordingly, the model chosen from the above-mentioned criteria was validated by Güner–Henry approach (decoy set) approach [60] for evaluating the robustness of the pharmacophore model on the basis of goodness of fit (GF) score in the range of 0 (null model) and 1 (ideal model) [58,61,62].

Given a particular dataset, pharmacophore validation is a quintessential criterion for

$$GF = \left(\frac{Ha}{4HtA}\right)(3A + Ht) \times \left\{1 - \frac{Ht - Ha}{D - A}\right\}$$

ferred to as active (A) mTOR inhibitors, with the remaining compounds being inactive. The Ligand Pharmacophore Mapping was executed for dataset screening, complemented with FAST algorithm and the GF score was calculated. The decoy set approach was instigated by evaluating the selected pharmacophore model on an external dataset (D) of 300 compounds obtained from the same biological assay [63]. This dataset was divided into 50 compounds exhibiting IC<sup>50</sup> < 100 nmol/L,

referred to as active (A) mTOR inhibitors, with the remaining compounds being inactive. The *Ligand Pharmacophore Mapping* was executed for dataset screening, complemented with *FAST* algorithm and the GF score was calculated.

#### *4.3. Virtual Screening of Marine Natural Product Library*

The validated pharmacophore model was escalated to screen the Marine Natural Product (MNP) library composed of 14,492 compounds (http://docking.umh.es/downloaddb, (accessed on 11 December 2020). The DS module *Ligand Pharmacophore Mapping* was employed in pursuit of identifying scaffolds mapping the pharmacophoric features. The compounds obtained from screening were further subjected to filtering by Lipinski's Rule of Five (Ro5) [64], Veber's rules [65] and pharmacokinetics by absorption, distribution, metabolism, excretion and toxicity (ADMET) for retrieval of drug-like compounds. Accordingly, the *Filter by Lipinski and Veber Rules* and *ADMET Descriptors* modules within DS were recruited for evaluation. The Ro5 and Veber's rules collectively oversee the physiochemical properties for efficient retrieval of compounds with molecular weight ≤500 kDa, number of hydrogen bond donors ≤5, compound's lipophilicity (log*P*) ≤5, number of hydrogen bond acceptors ≤10 and the number of rotatable bonds ≤10. The drug-like compounds so obtained were subjected to molecular docking with the mTOR kinase domain along with PP242 as reference inhibitor.

#### *4.4. Molecular Docking of Drug-Like Compounds with mTOR Kinase Domain*

Molecular docking methods explore the binding conformations adopted within the catalytic sites of macromolecular protein targets, thereby evaluating the vital phenomena for the intermolecular recognition process [66]. The drug-like compounds obtained from the above filtering criterion were subjected to molecular docking in Genetic Optimisation for Ligand Docking (GOLD) v5.2.2 automated docking software [67]. The compounds were evaluated on the basis of the GOLD default scoring function—Gold score [68]. This fitness score functions by scoring the summation of protein-ligand van der Waals interaction energy and hydrogen-bonding energy. The retrieved 3D crystallographic structure of mTOR (PDB ID: 4JT5) complexed with PP242 ATP-competitive inhibitor was prepared by utilizing the *Clean Protein* module in DS. Protein preparation was further carried out by adding the missing residues and hydrogen atoms. The water molecules were removed along with the bound PP242 ligand. Prior to docking, the performance of GOLD docking was assessed by re-docking the native PP242 ligand into mTOR. The drug-like compounds were subjected to minimization employing *Minimize Ligands* module in DS, preceding docking. Subsequently, molecular docking of drug-like compounds with mTOR was followed by applying the same docking parameters utilized for PP242 docking allowing for generation of 50 conformers per ligand. Clustering of the obtained conformations was carried out to obtain the largest cluster and each compound in the cluster was examined on the basis of higher Gold score than reference compound PP242, binding mode within the mTOR catalytic site and molecular interactions with the key residues of the mTOR kinase domain. The selected potential compounds acquired from this strategy were refined by MD simulations.

#### *4.5. Molecular Dynamics Simulation of Identified Hits*

MD simulation studies of compounds identified from above docking were extensively carried out to decipher the molecular dynamics in water and comprehend the interaction of hit compounds with vital residues of mTOR active site at the atomistic level. The docked structures of these marine hits in complex with mTOR were used as initial coordinates for simulations with GROMACS v2018 [69]. The mTOR and hit compounds were applied with CHARMm27 force field [70] and topologies generated with SwissParam [71] fast force field generation tool, respectively. The dodecahedron water box was utilized to solvate the systems with TIP3P water model and further neutralized with Na<sup>+</sup> counter ions. The steepest descent energy minimization was performed to dodge bad contacts and this was

followed by a two-fold equilibration. The NVT (constant number of particles, volume and temperature) equilibration at 300 K with a V-rescale thermostat supplemented with NPT (constant number of particles, pressure and temperature) equilibration at 1 bar pressure with a Parrinello-Rahman barostat [72] was orchestrated, each for 1000 ps. The LINear Constraint Solver (LINCS) [73] and SETTLE algorithms [74] were applied to monitor bond constrains and the geometry of water molecules. The long-range electrostatic interactions were calculated by means of Particle Mesh Ewald (PME) [75] and the equilibrated systems were subjected to production simulation runs of 30 ns. The acquired MD results were visualized and interpreted manually in visual molecular dynamics (VMD) [76] and DS. The binding free energy (BFE) scores were further computed for hit compounds by MM-PBSA executing *g\_mmpbsa* tool implemented in GROMACS [77]. For this purpose, 40 frames of mTOR-ligand complexes were selected evenly from the last 10 ns of MD trajectories and the BFE ∆*Gbind* was computed as per the below equation.

$$
\Delta G\_{bind} = \left( G\_{complex} - \left( G\_{protein} + G\_{ligand} \right) \right)
$$

#### **5. Conclusions**

A structure-based pharmacophore model, exploiting the crystal structure of mTOR serine/threonine kinase with its bound selective inhibitor Torkinib revealed fundamental pharmacophoric features required for mTOR inhibition at its ATP-binding pocket. A systematic virtual screening strategy with the model as a query, retrieved 3019 compounds from the Marine Natural Products library and subsequent filtering via Lipinski's Ro5, Veber's rule, and ADMET was able to acquire 135 drug-like compounds. Molecular docking of these compounds at the ATP-binding site of mTOR procured four marine compounds with higher dock scores than Torkinib and significant binding interactions with key residues of the pocket. These compounds also presented with good binding free energy scores and the energy contribution of essential residues unveiled that our hits can inhibit mTOR via Glu2190, Trp2239 and Cys2243 deemed requisite for selective inhibition over PI3K. The in silico ADME and toxicity analysis suggests three out of the four identified hits with acceptable pharmacokinetic profile for their in vivo disposition. Additionally, the biological origin of these hits was identified as marine fungus and sponge. Overall, we believe that our hits provide scaffolds for future drug optimization studies and, therefore, recommend these hit compounds from marine natural habitat as therapeutics for the treatment of cancer. The identification of marine-derived natural compounds embodies an essential platform for future drug discovery studies against various protein targets implicated in cancers.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/1424 -8247/14/3/282/s1, Figure S1: Overlay of the docked pose (orange) PP242 inhibitor of mTOR kinase with its crystal structure conformation (green) in PDB ID: 4JT5. Figure S2. (A) 3D and (B) representation of interaction between PP242 and catalytic residues of the ATP-binding pocket of mTOR. The compounds and interacting residues are represented as sticks. The hydrogen bonding interactions are shown as green dashed lines, the hydrophobic interactions are shown as pink and purple spheres and the van der Waals interactions are displayed as light green spheres. Table S1. The docking scores and intermolecular interactions of reference PP242 and Marine Natural Product (MNP) library compounds with mTOR kinase domain (PDB ID: 4JT5). Table S2. Binding free energy scores of identified Marine Natural Product (MNP) library hits with mTOR calculated through MM-PBSA methodology. Table S3. Lipinski's and ADME properties of identified marine hits to determine their drug-likeness. Table S4. Toxicity properties of identified marine hits to determine their drug likeness.

**Author Contributions:** Conceptualization, S.P.; methodology, S.P.; software, S.P.; validation, S.P.; formal analysis, S.P., V.K., S.R. and G.L.; investigation, S.P.; resources, K.W.L.; data curation, S.P., V.K. and K.W.L.; writing—original draft preparation, S.P.; writing—review and editing, S.P. and K.W.L.; visualization, S.P.; supervision, K.W.L.; project administration, K.W.L.; funding acquisition, J.C.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2020R1A6A1A03044344).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the article.

**Acknowledgments:** This research was supported by the Bio & Medical Technology Development Program of the National Research Foundation (NRF) and funded by the Korean government (MSIT) (No. NRF-2018M3A9A70-57263).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **A New Computer Model for Evaluating the Selective Binding Affinity of Phenylalkylamines to T-Type Ca 2+ Channels**

**You Lu <sup>1</sup> and Ming Li 2, \***


**Abstract:** To establish a computer model for evaluating the binding affinity of phenylalkylamines (PAAs) to T-type Ca 2+ channels (TCCs), we created new homology models for both TCCs and a L-type calcium channel (LCC). We found that PAAs have a high affinity for domains I and IV of TCCs and a low affinity for domains III and IV of the LCC. Therefore, they should be considered as favorable candidates for TCC blockers. The new homology models were validated with some commonly recognized TCC blockers that are well characterized. Additionally, examples of the TCC blockers created were also evaluated using these models.

**Keywords:** T-type calcium channel blocker; homology modeling; computer-aid drug design; virtual drug screening; L-type calcium channel

#### **1. Introduction**

As the only type of voltage-gated Ca 2+ channels that are activated at or near resting membrane potentials, T-type Ca 2+ channels (TCCs) play an important role in regulating [Ca 2+ ]<sup>i</sup> homeostasis in a variety of tissues, including pancreatic β-cells and tumor cells [1–4]. Therefore, TCC antagonists could be potentially useful for the treatment of chronic diseases associated with Ca 2+ dysregulation [5–7]. For this reason, it is imperative to develop more selective TCC antagonists for prospective clinical applications. Since many existing TCC blockers, such as mibefradil, also show inhibitory effects on L-type calcium channels (LCCs), the most important task in developing new TCC blockers is to enhance their selectivity to TCCs over LCCs. To achieve this, we established TCC– phenylalkylamine interaction models based on the specific amino acid sequences in the P-loop of TCCs, and α1C LCCs for characterizing the drug molecules' affinities for TCCs and LCCs, respectively.

TCCs have a close evolutional relationship with LCCs. A recent report from a cryoelectron microscopy study reveals that the frame of the α1G (Cav3.1) pore domain structure is similar to that of α1S (Cav1.1) [8]. This similarity allowed us to confidently adopt the global structure of the calcium channel CavAb model, constructed based upon Arcobacter butzleri crystallization [9], in the establishment of our TCC model. One of the most remarkable differences between all types of TCCs and LCCs is a lysine residue located adjacent to the critical glutamic acid/aspartic acid residue in domain III. The existence of a positively charged lysine (K 3p49 ) may swing the aspartic acid (D3p50 ) away from the center of the calcium filter and change the preferred calcium ion and drug binding sites from domains III and IV for LCCs to domains I and IV for TCCs. Therefore, we used the ZMM molecule modeling program [10–12] to create four-domain TCC models, in which the binding affinities of drugs to TCCs and α1C LCC were determined by scoring their free energy in binding to the channels [13].

The new TCC models are adopted from a drug–protein interaction framework for modeling CavAb blocking by phenylalkylamines (PAAs) [9]. This is rational because many TCC blockers are PAAs or their derivatives, and because PAAs block CavAb [9,14].

**Citation:** Lu, Y.; Li, M. A New Computer Model for Evaluating the Selective Binding Affinity of Phenylalkylamines to T-Type Ca 2+ Channels. *Pharmaceuticals* **2021**, *14*, 141. https://doi.org/10.3390/ ph14020141

Academic Editor: Osvaldo Andrade Santos-Filho Received: 27 January 2021 Accepted: 8 February 2021 Published: 10 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

It is proposed that PAAs bind to LCCs in an inverse V-shaped configuration, with the ammonium group towards the P-helices, and the nitrile group bound to the calcium ion coordinated by the selectivity filter glutamates in domains III and IV of the LCC [15]. We reason that this is also true for TCC blockers, except that the calcium ion is coordinated in the cavity between domains I and IV, since the depolarization confirmed that the pore domains of CaV3.1 and CaV1.1 are superimposed [8]. As a result, the two rings of the flexible PAA molecule [15] will make hydrogen bonding contacts with the mobile side chains of relevant amino acids from domains I and IV of TCCs. This strategy allowed us to create computer models for simulating the interactions between drugs and channel receptors for LCC and TCCs, respectively.

#### **2. Results**

#### *2.1. Homology Modeling of TCCs and α1C LCC*

Using the bacterial calcium CavAb open channel 3D structure as the input, ZMM generated the first template of the calcium channel, which was then modified with S5-Ploop-S6 segments of α1C, and α1G, α1H, and α1I (Table 1) to create corresponding protein structures of Cav1.2 LCC and Cav3.1, Cav3.2, and Cav3.3 TCCs, respectively (Figure 1A–D). For cross-validation of ZMM-generated 3D structure models, we also performed ab initio modeling of α1C and α1G calcium channels. Since there is a considerable overlap of PAA inhibition between LCC and TCCs [9], the allosteric structures of LCC and TCCs are more likely to be similar. Comparing two different homology modeling tools, ZMM generates more consistent 3D models of the domain III S5-P-loop-S6 segment of α1C and α1G (Figure 1E,F) than the ab initio method (Figure 1G,H). *α* α α α α α α α α

α α α α α α α α α α α α α α **Figure 1.** Top views of homology modeling results from ZMM for α1C, α1G, α1H, and α<sup>1</sup> I. (**A**–**D**) The structures of α1C, α1G, α1H, and α<sup>1</sup> I, respectively; blue, green, brown, and yellow represent domains I, II, III, and IV, respectively; the four selectivity-determining amino acids (glutamic acid or aspartic acid) in the P-loop are colored red and displayed as spheres; ZMM generates more consistent 3D structure than the ab initio modeling method for α1C and α1G; (**E**) the predicted 3D structure of the α1C domain III generated by ZMM, the glutamic acid is represented by red spheres; (**F**) the predicted 3D structure of the α1G domain III generated by ZMM, the lysine is represented by red spheres; (**G**) the most representative structure selected by Calibur clustering analysis [16] of α1C domain III, the glutamic acid is represented by red spheres; (**H**) the most representative structure selected by Calibur clustering analysis of α1G domain III, the lysine is represented by red spheres.

#### *2.2. Further P-Loop Remodeling of TCCs*

After determining the globe structure of TCC 3D models, we focused on the variability of the P-loop structure, which is the major drug–ligand interaction segment. The Rosetta P-loop remodeling module [17] was utilized to estimate the variability of P-loop 3D structures on every domain of TCCs. After inputting a perturbation to the original structure, the remodeling process was conducted by sampling the possible locations of a given length of an amino acid sequence in three-dimensional space. Using the ZMM generated α1G structure as the reference, the energy-based clustering method [18] was used to determine the P-loop remodeling results with the lowest root-mean-square-displacement (RMSD) score. We found that for α1G, domain II had a clear variation between two different homology modeling methods in sample sizes 500 and 20,000 (Figure 2). It showed that the central P-loop helix segment is in the horizontal position rather than the diagonal position found in other domains. As a result, the selectivity-determining glutamic acids E2p50 may have a larger vertical distance from other glutamic acids/aspartic acids (E1p50, D3p50, and D4p50) in α1G TCCs. This may exclude glutamic acid E2p50 as a Ca2+ binding candidate, leaving the Ca2+ to bind either E1p50 to D4p50 or D3p50 to D4p50 in TCCs. Additionally, to validate the normality of the remodeling data, we conducted a nonparametric test for the α1G P-loop remodeling data and confirmed that all the sampling processes (500 and 20,000) came from the same distribution (see Supplementary Table S2, Supplementary Figure S2 for statistical results).


**Table 1.** Comparison of numerical results of P-loop electrostatic potential at the four different domains with different lengths of amino acid sequences. TCC: T-type calcium channel.

### *2.3. Local Electrostatic Potentials of the Selective P-Loop of TCC Domains and the Impact of K3p49*

A previous study indicated that when a calcium ion enters the selectivity filter region of a LCC, it binds to the selectivity-determining glutamic acids (E3p50, E4p50) in domains III and IV [15]. Consequently, the phenylalkylamine molecules will bind to domains III and IV due to the interaction between the nitrile nitrogen and Ca2+ [15]. In contrast, all TCCs have a lysine (K3p49) located at the 5′ end adjacent to D3p50 in domain III. It is reported that the replacement of lysine (K3p49) with Phe or Gly causes the activation curve to shift to the right [8], which indicates that the lysine in the position adjacent to aspartic acid (D3p50) plays a significant role in the kinetic/dynamic mechanism of Ca2+ interaction with the inner environment of the central cavity of TCCs. This positively charged lysine alters the negative charge field distribution of D3p50 to attract Ca2+ (Supplementary Figure S1). It is possible that the lysine (K3p49) swings aspartic acid (D3p50) away from the original Ca2+ binding position, thus causing the Ca2+ to bind glutamic acid or aspartic acid in other domains, probably to domains I and IV since domain II has a configuration deviation. As a result, the phenylalkylamine may also switch its binding region from domains III and IV to domains I and IV.

α

α

α

α

α α α **Figure 2.** Comparison of the P-loop conformation differences before and after Rosetta P-loop remodeling. (**A**) Homology modeling of P-loop structures of α1G domain I (**A1**), domain II (**A2**), domain III (**A3**), and domain IV (**A4**) generated by ZMM; (**B**) Rosetta P-loop remodeling results of α1G domain I (**B1**), domain II (**B2**), domain III (**B3**), and domain IV (**B4**) with the sampling size equal to 500; (**C**) Rosetta P-loop remodeling results of α1G domain I (**C1**), domain II (**C2**), domain III (**C3**), and domain IV (**C4**) with the sampling size equal to 20,000.

To determine the effect of lysine (K3p49) on overall electrostatic potential (E) for given TCC homology models, we calculated the electrostatic potential (E) for the tailed P-loop of domain I to IV. Table 1 shows that the combined electrostatic potential (Ecoul) becomes more negative in each domain as the number of testing amino acids is reduced from seven to five. In TCCs, Ecoul for domain I is more negative than that for domain IV in the five amino acid-reduced sequence, indicating a possible switching of the Ca2+ binding site from domains III and IV to domains I and IV.

We used Coulomb's electric force equation to quantitively analyze the influence of lysine on the electrical attraction force between Ca2+ and aspartic acid (D3p50). According to the equation in Section 4.2, lysine (K3p49) has the least effect on the Ca2+-D3p50 attraction when K3p49 is located on the opposite side of the Ca2+ and when D3p50 is at the center. When *a* = 4.3 Å and *b* = 3.8 Å [8] (calculated in PyMOL.2.3.3 for ZMM results), Coulomb's force equation (found in Section 4) yields: F(Ca, D, attraction) = 2.489 × 10−<sup>9</sup> N and F(Ca, K, repellent) = 0.702 × 10−<sup>9</sup> N; thus, lysine, at a minimum, reduces the attraction force between Ca2+ and aspartic acid by more than 28%. The attraction force between Ca2+ and D3p50 will reduce further or reverse into a repellent force as the distance from the Ca2+ to K3p49 decreases; therefore, the preferred binding position of Ca2+ will likely be switched to domains I and IV in TCCs. This limits the PAA binding region on TCCs to domains I and IV. We could use the amino acid structure of domains I and IV to evaluate the affinities of the PAAs (and their derivatives) for TCCs (using models established for α1G, α1H, and α1I) and use the amino acid structure of domains III and IV for evaluating their binding affinities to LCC (using the model of α1C).

### *2.4. Model Predictions and Vina Screening Output of Some Current T-Type Ca2+ Channel Blockers*

α

α α α

Mibefradil is reported to have an inhibitory effect on both LCC and TCCs [20]. The α1C model predicts that one hydrogen atom from the nitrogen (N3) on the cyclopentadiene connects to methionine (M4i27) on domain IV of α1C LCC, as shown in Figure 3A,B. This is consistent with the prediction of another model in a previous study [9]. In contrast, the α1C homology model does not predict that hydrogen bonds to NNC 55-0396. NNC 55-0395 inhibits both L- and T-type calcium channels [14]. The α1C homology model predicts that NNC 55-0395 has one hydrogen bond that connects nitrogen (N3) to the glycine (G3p49) on P-loop domain IV. For NNC 55-0397, the α1C homology model also predicts that RO 40-5966, a hydrolyzed metabolite of mibefradil [20], has one hydrogen atom from the hydroxy group of the benzene ring bound to the glycine (G4p49) at domain IV. No hydrogen bond has been found between the LCC and the TCC blocker SKF-96365. α α α α α

α α α α **Figure 3.** Models of ligand–receptor interactions of mibefradil and NNC 55-0396. (**A**) The predicted binding sites of mibefradil on α1C. The H-bond formed between the ammonia (N<sup>3</sup> ) on the cyclopentadiene of mibefradil and methionine (M4i27) on domain IV. The relative locations of surrounding amino acid residues of the α1C L-type calcium channel (LCC) are shown by the arch–dash symbols. (**B**) The predicted 3D binding sites of mibefradil on α1C domain IV from a side view. Red spheres represent the position of glutamic acid E4p50. Mibefradil is represented by the green ring structure. (**C**) The predicted binding sites of NNC 55-0396 on α1G domain IV. The H-bond formed between the central ammonium of NNC 55-0396 and asparagine (N1i20) on domain I. (**D**) The predicted 3D binding sites of NNC 55-0396 on α1G. NNC 55-0396 is represented by the green ring structure. Red spheres represent the position of glutamic acid E4p50. For A and C, carbon, nitrogen, oxygen, and fluorine elements are represented by black, blue, red, and green, respectively; for B and D, the blue and orange ribbon helices represent S5 and S6, respectively. The ribbon helix structures linking S5 and S6 are P-loops. The yellow ball represents the position of the calcium ion.

Using TCCs as templates, we have revealed some current TCC blockers of α1G, α1H, and α1I. The α1G model predicts that the fluorine atom from the compound NNC 55-0395 forms a halogen bond to glycine (G1p51) in domain I. Our α1G model also predicts a binding site of NNC 55-0396 to asparagine (Nli20) in domain I (Figure 3C,D). For NNC 55-0397, the

fluorine atom from the compound forms a halogen bond with valine (V1p46) at the P-loop of domain I. For mibefradil, one oxygen atom from the side chain of the compound forms hydrogen bonds with asparagine (N1o4) at S5 of domain I. For RO 40-5966, one hydrogen atom from the nitrogen (N3) on the cyclopentadiene forms a hydrogen bond with alanine (A1i27) in domain I. For SKF-96365, the center oxygen atom forms a hydrogen bond with asparagine (N1o4) at S5 of domain I.

Our α1H model predicts that one hydrogen atom from the nitrogen (N3) on the cyclopentadiene of NNC 55-0395 interacts with valine (V1p46) at the P-loop of domain I to form a bond. The fluorine atom from NNC 55-0396 forms a halogen bond with isoleucine (I1i8) from the α1H S6 of domain I. For NNC 55-0397, one hydrogen atom from the nitrogen (N3) on the cyclopentadiene forms a bond to asparagine (N4p51) at the P-loop of domain IV. The hydrogen atom from the nitrogen (N3) on the cyclopentadiene of mibefradil finds asparagine (N1o4) to form a bond at S5 of domain I. For RO 40-5966, the fluorine atom from the compound forms a halogen bond with histidine (H4i29) at S6 of domain IV. One oxygen atom from the side chain of SKF-96365 forms a hydrogen bond with asparagine (N1o4) at S5 of domain I.

Our α1I homology model predicts that NNC 55-0395 forms a halogen bond between the fluorine atom from the compound and isoleucine (I1i8) at S6 of domain I. The hydrogen atom from NNC 55-0396 forms a bond to asparagine (N1o4) at S5 of domain I. For NNC 55-0397, one hydrogen atom from the nitrogen (N3) on the cyclopentadiene interacts with valine (V1p46) at the P-loop of domain I to form a bond. For mibefradil, there is a halogen bond formed between a fluorine atom from the compound and a hydrogen atom from asparagine (N4p53) at the P-loop of domain IV. For RO 40-5966, one hydrogen atom from the ammonia on the cyclopentadiene interacts with asparagine (N1o4) at S5 of domain I to form a bond. Our model does not predict the hydrogen bond formed when docking SKF-96365 to α1I.

A comparison of the predicted binding affinity K<sup>d</sup> and experimental measurements of IC<sup>50</sup> for given TCC blockers are listed in Supplementary Table S4. Table 2 summarizes the predicted binding affinity results for all existing TCC blockers.

#### *2.5. Evaluation of New Compounds*

The Vina [13] models were employed for evaluating the binding affinity of the testing compounds. We randomly selected 300,000 compounds from PubChem and used these as the database to train our recurrent neural networks (RNNs) [21] with the given compound properties.

After performing virtual screening, we found that the compounds TC 7, TC 4, and TC 2 satisfied our screening criteria for α1G, α1H, and α1I, respectively. Compound TC 7 has the highest binding affinity, as well as a lower (water–octanol partition coefficient) logP and a higher Quantitative Estimation of Drug-likeness (QED) than existing TCC blockers. The 3D binding plots between TC 7 and α1G are shown in Supplementary Figure S3. The predicted binding affinities between existing TCC blockers and screened compounds on TCCs and α1C LCC are shown in Figure 4. Our results show that these screened compounds have smaller logP and Synthetic Accessibility Scores (SAS) but larger QED values than those of selected TCC blockers (as seen in Table 3). More structures and chemical properties for the 13 identified compounds can be found in Supplementary Figures S5–S17.


**Table 2.** Predicted Gibbs free energy ∆G of phenylalkylamines (PAAs) on calcium channels (N/A: not available).

**Figure 4.** Predicted Gibbs free energy of select T-type Ca2+ channel (TCC) blockers and computerdesigned compounds of different receptors. Blue—α1C, brown—α1G, gray—α1H, and yellow—α<sup>1</sup> I.


**Table 3.** The chemical properties of computer-designed compounds and selected TCC blockers. SAS: Synthetic Accessibility Scores, QED: Quantitative Estimation of Drug-likeness.

#### **3. Discussion**

A TCC (Cav3.1) 3D structure has already been modeled with cryo-electron microscopy [8]; however, this structure is constructed based on a splice variant containing a deletion of 133 amino acids within the I-II linker. Electrophysiological characterization of these variants (Cav3.1-∆8b) shows 1.5-2-fold conductance increases when compared with the full-length form in human and rat preparations. Both activation and steady-state inactivation curves are shifted in the human preparation [8]. In addition, the pore diameter estimated from Cav3.1-∆8b is smaller than the biophysical measurement [8]. These alterations in TCC electrophysiological properties suggest that the conformation of the cryo-electron microscopy structure is not the same as the full-length Cav3.1 TCC. Therefore, the 3D structure of Cav3.1-∆8b may not be the most suitable template for the general modeling of TCCs, especially for PAA binding, which is highly dependent on the position of Ca2+ interacting with the selectivity filter of TCCs. In this study, we chose to use CavAb as the model template since this channel is blocked by PAA and therefore is suitable for establishing a model for evaluating PAAs that inhibit TCCs selectively over LCCs.

Increasing evidence indicates the pathological role of TCCs in the progression of different diseases [6]. It is crucial to develop selective TCC blockers to establish new treatments for these diseases. Unfortunately, lacking the TCC X-ray crystallization structure hampers the progress of creating new TCC blockers. In practice, it is very difficult to find or design a compound that selectively blocks TCCs but not LCCs, since most current TCC blockers exhibit a certain level of inhibitory effects on LCCs. For example, mibefradil, the first launched TCC inhibitor, was quickly recognized to cross inhibit LCC [20]. Here, we provide a new strategy by which the specificity of candidate compounds for binding TCCs but not LCCs can be pre-screened with new computer-based models. This is desirable for designing and developing compounds that more selectively block TCCs than LCCs.

We chose to build models for the interaction between TCCs and PAAs since the binding mechanism of these compounds to LCCs has been studied extensively [9,15,22–26]. Based upon the critical single amino acid lysine (K3p49) difference between LCCs and TCCs,

we have created a strategy that can distinguish the affinity of PAAs to TCCs and LCCs, respectively. The models have been validated by measuring the affinities of existing TCC blockers to LCCs and TCCs. We also used these models for evaluating the specificities of novel PAAs and phenylalkylamine derivatives in terms of their binding affinities to TCCs and LCCs.

Using ZMM, we simultaneously created α1G, α1H, and α1I TCC and a1C LCC structures with four domains, each containing three segments: segment 5, a P-loop, and segment 6. In contrast, the ab initio modeling method failed to produce a suitable calcium channel structure compared to ZMM.

The selectivity of the calcium channel is dependent on the critical glutamate residues located in the selectivity filter of the P-loop of the α<sup>1</sup> subunit in each domain. In this region, there are two negatively charged glutamic acid residues likely to attract one Ca2+ in the space close to domains III and IV [15]. When a phenylalkylamine molecule approaches a calcium channel from the cytoplasmic side, its nucleophilic nitrile nitrogen reaches the Ca2+, while the other parts of the molecule form affiliated interactions with the amino acids in the P-loop and segments 5 and 6 in domains III and IV of the calcium channel. This causes a physical blockage of ion flow through the channel. In the case of TCCs, there is a lysine (K3p49) located adjacent to D3p50 in domain III, and the ionized electric potential distribution of the aspartic acid is altered by lysine, which attenuates the electric attraction of aspartic acid (K3p49) to Ca2+ at the minimum binding distance (4.3 angstroms) and may swing K3p49 away from the selectivity filter. Based on this analysis, we suggest that Ca2+ will not bind to domain III but to domain I of TCCs. This prediction is consistent with the Cav3.1 structure estimated with cryo-electron microscopy [8], which showed the electron density of the top Ca2+ ion is closest to Glu354 of Cav3.1 (E1p50, Table 4). Additionally, Rosetta P-loop remodeling shows that the P-loop of domain II is in a more horizontal confirmation than that of other domains, rendering the glutamic acid (E2p50) further away from the Ca2+ binding site. Since the movements of PAAs will follow the location of Ca2+ docking, our models are built for evaluating the affinity of candidate compounds binding to domains I and IV. The compounds that are predicted to have a higher affinity to bind domains I and IV of TCCs but not domains III–IV of LCC (α1C) are considered to be ideal selective TCC blocker candidates. This strategy screens out the compounds that are unlikely to bind LCC and TCCs as well as the compounds that are likely to bind both LCC and TCCs. To test whether Ca2+ docking is consistent on domains I and II across different TCCs, a molecular dynamics study should be conducted with modified membrane conditions and simulation environments [27].



734

V2i8, L3o7, N3p54, S3i15, Y4o16, K4o29, Q4p45, F4i1, V4i2, I4i19, V4i23. a,b

 are labeled according

 to the alignment

 of the outer helix, the P-loop, and the inner helix of the KcsA structure

 [15].

Residue sequences

Previous studies suggest that the nitrile and isopropyl groups in devapamil and some other PAAs serve to guide the drug to the position of Ca2+; this function persists if the nitrile is replaced with other elements with high electronegative potentials, such as oxygen or sulfur [15]. In many molecules discussed here, including mibefradil, the nitrile is replaced with a methoxy acetyl side chain or a similar side chain with a high electronegative potential. These molecules behave presumably like those of molecules with nitrile in their alkaline chain. Some of the molecules, such as RO 40-5966 and SKF-96365, do not share the binding mechanism described by our model, and therefore their inhibitory effects on Ca2+ channels may not be explained by our new models. For example, by using the input template SKF-96365, we obtained 14 unique structures (as seen in Supplementary Table S3, Figure S4), which had negative binding affinities to our TCC models.

Although our models are designed to select compounds that are likely to bind domains I and IV of TCCs, this does not exclude the possibility that PAAs or their derivatives might inhibit TCCs via binding to domains I and II, domains II and III, or even domains III and IV. The goal of our models was to increase the likelihood of success in screening selective TCC blockers based on their chemical structures.

Our α1C model has a similar channel pore size (selectivity filter region) as CavAb [9]; however, α1G, α1H, and α1I may have smaller diameters than CavAb, since the unitary conductance of TCC currents is smaller than that of LCC. Further statistical analyses of P-loop remodeling data show that a minimal structural difference exists in the P-loop region remodeling data (see Supplementary Materials for details of the statistical analysis).

The Vina screening results of α1C identify no binding location for NNC 55-0396 or SKF-96365. The predicted α1C binding amino acid for NNC 55-0397, as well as mibefradil, matches the experimental results, which show the inhibitory effect of NNC 55-0397 and mibefradil on LCCs. Although RO 40-5966 has a lower ∆G than mibefradil when binding to α1C, the predicted binding location is closer to the center of the channel filter region than mibefradil, which indicates a stronger blocking effect on the rate of Ca2+ influx than mibefradil.

Although the K<sup>d</sup> values predicted by Vina have some gaps compared to the experiment data, they do follow the same order of magnitude (Supplementary Table S4). To obtain a more accurate Gibbs free energy for PAAs binding to TCCs, at least two consecutive steps must be conducted: first, the flexible docking process [28]; second, the free energy calculation between ligand and receptor [29]. These two steps require an extensive computational cost and the final K<sup>d</sup> value is very sensitive to the initial input of the receptor structures. Recently deposited human α1G structures offer a good template for developing TCC blockers [8]; however, they have some uncommon regions missing, which could affect PAA binding. Therefore, we argue that it is less likely that the Gibbs free energy of mibefradil/NNC 55-00396 between the prediction and the experiment is matched by choosing different docking programs or conducting a molecular dynamics simulation to find the free energy.

Our work only focuses on the first step of the drug development process in silico, providing a strategy for predicting the comparative potency of candidate compounds to TCCs versus LCCs. Neither are used for evaluating the pharmacological effects of these compounds on other types of cation channels. Since the strength of the pharmacological effects of PAAs and their derivatives on blocking calcium channels are increased by the appearance of a calcium cation in the channel pore [15], it is unlikely that these PAAs and their derivatives will exhibit a strong inhibitory effect on other cation channels.

#### **4. Materials and Methods**

#### *4.1. Homology Modeling of the α<sup>1</sup> Subunit*

Three classes of calcium channel families have been discovered: CaV1.X, CaV2.X, and CaV3.X. The X represents the subdivisions of the sequence homology of the α<sup>1</sup> subunit in each class. The models of drug–channel interactions built upon the structural differences in the relevant S5, P-loop, and S6 regions of α1G, α1H, α1I TCCs, and α1C LCC, respectively. The protein templates were obtained from BAM [30] using truncated inputs of human α1G, α1H, α1I, and α1C (UniProt id: O43497) amino acid sequences (see Table 4 for details). The crystallization structure (PDB id: 5kmh) for the depolarization status of the calcium channel protein, originally extracted from Arcobacter butzleri [9], was employed as the structure template of our model. The multi-domain protein structures of human α1G, α1H, α1I, and α1C were built using the ZMM molecular modeling software. The forcefield of specific amino acids was simulated by using the Assisted Model Building with Energy Refinement (AMBER) program. The final structure of the target peptides was optimized by using the Monte Carlo minimization protocol. The maximum iteration time for finding the global minimum was set to 5000. During the energy optimization, structural similarity between target and template was maintained by a flat-bottom parabolic energy penalty function that allows for penalty-free deviations of alpha-carbons up to 1 atom distance from their respective positions in the template, and a penalty was imposed with a force constant of 10 kcal mol-1A-2 for larger deviations [15]. The homology models for human α1H (UniProt id: O95180), α1I (UniProt id: Q9P0X4), and α1C (UniProt id: Q13936) were also built with this method.

#### *4.2. Local Electrostatic Potential Calculation*

To calculate the electric double layer-related local electrostatic potential while a channel protein interacted with a surrounding water molecule, we generated the corresponding meshes using MSMS (v2.6.1) [31] and set the probe radius to 1.4 and the density to 3.0 for quality control. For truncated amino acid sequences, the local electrostatic potential/binding energy is derived from the summation of the solvation energy and Coulomb energy, i.e., Gcomplex = Gsolution + GCoulomb. The GCoulomb for LCC and TCCs in the P-loop region was calculated by using PyGBe with pre-defined parameters (Supplementary Table S1). To compare the influence provided by a single lysine, we used Coulomb's law to calculate the electric attracting force between Ca2+ and aspartic acid:

$$\mathbf{F}\_{\text{(Ca,D)}} = \mathbf{k}\_{\text{e}} \frac{\mathbf{q}\_{\text{Ca}} \mathbf{q}\_{\text{D}}}{\mathbf{a}^2}$$

where k<sup>e</sup> is Coulomb's constant 8.99 × 109 N·m<sup>2</sup> ·C −2 ; a is the distance between Ca2+ and aspartic acid; q is the point charge for Ca2+, aspartic acid (D), and lysine (K), respectively. The effect of the repellent force on Ca2+ by lysine in the direction of the Ca2+ and the aspartic acid attracting force is defined by

$$\mathbf{F}\_{\text{(Ca,K)}} = \mathbf{k}\_{\text{e}} \frac{\mathbf{q}\_{\text{Ca}} \mathbf{q}\_{\text{K}}}{\mathbf{r}^2} \cdot \cos \theta\_{\text{e}}$$

where θ is the angle between the lines from Ca2+ to aspartic acid and from Ca2+ to lysine; r is the distance between Ca2+ and lysine, which is calculated by

$$\mathbf{r} = \mathbf{a} \cdot \cos \theta \pm \sqrt{\mathbf{b}^2 - (\mathbf{a} \cdot \sin \theta)^2} \dots$$

where b is the distance from aspartic acid to lysine. When the angle between Ca2+ and aspartic acid, and lysine and aspartic acid (φ) is less than 90◦ , r = a· cos θ + q b <sup>2</sup> − (a· sin θ) 2 ; when φ is equal to 90◦ , r = a· cos θ; when φ is larger than 90◦ , r = a· cos θ – q b <sup>2</sup> − (a· sin θ) 2 . Therefore, the electric force between Ca2+ and aspartic acid is

$$\mathbf{F}\_{\text{(Ca, D, K)}} = \mathbf{F}\_{\text{(Ca, D)}} - \mathbf{F}\_{\text{(Ca, K)}} = \begin{vmatrix} \mathbf{q}\_{\text{Ca}} \mathbf{q}\_{\text{D}} \\ \mathbf{a}^2 \end{vmatrix} - \mathbf{k}\_{\text{e}} \frac{\mathbf{q}\_{\text{Ca}} \mathbf{q}\_{\text{K}}}{\mathbf{r}^2} \cdot \cos \theta.$$

#### *4.3. Ab Initio Modeling*

The ab initio modeling modules from Rosetta were employed to find the threedimensional structure of target fragmental peptides by sampling and assembling a large

candidate pool containing 22,000–27,000 decoy structures for every inputted amino acid sequence [32]. The output results of ab initio modeling were analyzed using the Calibur and energy-based clustering methods.

#### *4.4. P-Loop Remodeling*

The P-loop region of α1G was remodeled using a Rosetta loop modeling module combined with the FastRelax protocol. Twenty-seven amino acids in the P-loop were selected from domains I, II, III, and IV of the TCCs. The modeling used phenylalanine as the starting amino acid and tyrosine, tryptophan, proline, and isoleucine as the ending amino acids. The effective sample size used for subsequent statistical analysis was validated by two groups of data for every remodeled domain. The first group contained 500 output structures, and the second group had 20,000 output structures. The results were analyzed using the Calibur and energy-based clustering methods.

#### *4.5. Compound Generation*

We used the de novo drug generation package "chemical vae" developed by Gomez-Bombarelli et al. [21] to create a data-driven RNN for new compound production. The dataset we used to train the RNN was prepared by randomly sampling approximately 250,000 compounds from PubChem. The maximum length of encoding for the SMILES-based compounds was set to 120 characters. To analyze the compounds using the RNN, one fully connected layer of width 200 was used. To convert a predicted compound back to the original data type, three layers of gated recurrent units with a hidden dimension of 500 were used. The variational loss of the RNN was annealed according to a sigmoid schedule after 35 epochs, running for 130 epochs while property prediction training the RNN, such that the RNN trained on the PubChem data set with objective properties including: logP, SAS [33], and QED [34]. We kept the other hyperparameters to train the RNN unchanged from the reference [21]. To transfer the predicted compound back to SMILES-based data, we set the Gaussian noise value to 5 and the iteration time to 1000. Once the 2D structure had been obtained, we converted it into a 3D structure via the online program Frog 2.1 [35,36]. The program OpenBabel 2.4.1 [37] was used to add the hydrogen atom and set the pH equal to 7.35 for select compounds.

The 2D structure of mibefradil was employed as a redesigned template for new compounds. Based on its structure, we recreated 129 PAAs and their derivatives. Their corresponding 3D structures (involving up to 800 isomers) were created by Frog 2.1. We used the same program to find the 3D structures for NNC 55-0395, NNC 55-0396, and NNC 55-0397, and combined them with SKF96365 and RO 40-5966, whose 3D structures were downloaded from PubChem, for use as reference compounds for testing and validating the faithfulness of our TCC models.

Some of the candidates of screening compounds may contain oxygen, which replaces the role of nitrile; this structural formula has been reported in certain PAAs such as falipamil, BRL-32872, and tiapamil [15].

#### *4.6. Virtual Drug Screening*

Virtual drug screening was conducted by using AutoDock Vina [13] with user-defined configuration scripts on the Tulane supercomputer Cypress. The search box was placed in the center of the protein model. The number of mesh elements in the X, Y, and Z directions was set to 60, 124, and 102 for α1C and 58, 48, and 50 for α1G, α1H, and α1I, respectively, when simulating the ligand–receptor interaction for existing TCC blockers. The number of mesh elements in the X, Y, and Z directions was increased to 126 when conducting virtual screening for newly designed compounds. To achieve repeatable docking results using Vina, the seed number was fixed at –1460306363. As the grid number for every direction was set to the maximum, Vina had to search a very large three-dimensional space. To find the local minimum, the exhaustiveness was set to 2000 for new compound screening cases and 8 for existing TCC blocker screening cases. The number of predictable binding models

expressed as the output was limited to 3 for new compound screening cases and 20 for existing TCC cases.

The Vina output results were checked using PyMOL to ensure the binding locations for existing and newly designed compounds. The predicted binding affinity for the testing compound was calculated as:

$$\mathbf{K\_d} = \exp\left(\frac{-\Delta\mathbf{G} \cdot \mathbf{kcal} \cdot \mathrm{mol}^{-1}}{0.001986 \cdot \mathrm{kcal} \cdot \mathrm{mol}^{-1} \cdot \mathrm{K}^{-1} \cdot 310 \,\mathrm{K}}\right) \,\mathrm{s}$$

where ∆G is the Gibbs free energy predicted by Vina.

The 2D ligand–receptor interaction plot was created using LigPlot<sup>+</sup> [38].

#### *4.7. Data Analysis*

The Anderson–Darling normality test and the Kruskal–Wallis one-way ANOVA test (Supplementary Table S2) were conducted on the generated homology modeling data from Rosetta ab initio modeling and P-loop remodeling in Anaconda Spyder (3.2.8) using a Python 3.6 environment.

#### **5. Patents**

All the new identified compounds in this study are patented by the Office of Technology Transfer and Intellectual Property Development at Tulane University (Patent ID: US62/859,519).

**Supplementary Materials:** The following are available online at https://www.mdpi.com/1424-824 7/14/2/141/s1. Figure S1. The negatively charged lysine affects the electric potential distribution of aspartic acid in the x-y plane. (A) The Dist. of E<sup>D</sup> (0,0) without K; (B) the Dist. of E<sup>D</sup> (0,0) after adding the lysine E<sup>K</sup> (3.8, 0), Figure S2. The predicted versus theoretical RMSD plots for group sampling size 500 and 20,000. (A,B,C,D) Domains I to IV for 500 sampling size group; (E,F,G,H) Domains I to IV for 20,000 sampling size group Figure S3. The predicted 3D binding plots between TC 7 and α1G. An alkyl bond (4.31 angstroms) has been formed between TC 7 (green) and the sidechain of V1i24 (red) at domain I. The sidechain of E1p50 is colored as blue and one Ca2+ is colored as yellow; Figure S4. Based on the 2D structure of SKF-96365, 10 structures have been found using the Deep-Learning based de novo drug design approach, Figure S5. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 1, Figure S6. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 2, Figure S7. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 3, Figure S8. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 4, Figure S9. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 5, Figure S10. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 6, Figure S11. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 7, Figure S12. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 8, Figure S13. 2D structure and chemical properties of redesigned phenylalkylamine analog, TC 10, Figure S14. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 11, Figure S15. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 12, Figure S16. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 13, Figure S17. 2D structure and chemical properties of the redesigned phenylalkylamine analog, TC 15. Table S1. Numerical parameter settings for running PyGbe, Table S2. Normality test for P-loop remodeling data (α1G) from two groups with different sampling sizes, Table S3. The structures and properties of computer-designed compounds using Deep-Learning (D: distance; C: count; F: frequency). Table S4. Predicted binding affinity (Kd) by Vina versus experimental measurement of IC<sup>50</sup> of given TCC blockers (unit: micromolar).

**Author Contributions:** Conceptualization, Y.L. and M.L.; methodology, Y.L. and M.L.; software, Y.L.; validation, Y.L. and M.L.; formal analysis, Y.L.; investigation, Y.L.; resources, Y.L.; data curation, Y.L.; writing—original draft preparation, Y.L.; writing—review and editing, Y.L. and M.L.; visualization, Y.L.; supervision, M.L.; project administration, M.L.; funding acquisition, M.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported, in part, through the use of high-performance computing (HPC) resources and services provided by Technology Services at Tulane University, New Orleans, LA.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** Y.L. thanks the interdisciplinary Ph.D. Program in Aging Studies at Tulane University for their support of this research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Pharmaceuticals* Editorial Office E-mail: pharmaceuticals@mdpi.com www.mdpi.com/journal/pharmaceuticals

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-5383-2