Molecular modeling technology studies of novel pyrazoylethylbenzamide derivatives as selective orexin receptor 1 antagonists

Molecular modeling technology studies of novel pyrazoylethylbenzamide derivatives as selective orexin receptor 1 antagonists

Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17 Contents lists available at ScienceDirect Journal of the Taiwan Institute of C...

6MB Sizes 0 Downloads 70 Views

Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

Contents lists available at ScienceDirect

Journal of the Taiwan Institute of Chemical Engineers journal homepage: www.elsevier.com/locate/jtice

Molecular modeling technology studies of novel pyrazoylethylbenzamide derivatives as selective orexin receptor 1 antagonists Qunlin Li, Cuihua Zhang, Yujie Ren∗ College of Chemical and Environmental Engineering, Shanghai Institute of Technology, Shanghai 201418, China

a r t i c l e

i n f o

Article history: Received 14 January 2019 Revised 26 March 2019 Accepted 28 March 2019 Available online 10 May 2019 Keywords: Orexin receptor (OXR) Insomnia Pyrazoylethylbenzamide derivatives Computer-aided drug design Virtual screening

a b s t r a c t Orexin neuropeptide signal is the core promoter of arousal, and orexin signaling is activated by orexin receptor 1 (OX1 R), which makes OX1 R an attractive target for the treatment of insomnia. In this study, OX1 R quantitative structure-activity relationship (QSAR) models (CoMFA: q2 = 0.788, r2 = 0.988; CoMSIA: q2 = 0.783, r2 = 0.953) were constructed based on the pIC50 values of 36 pyrazoylethylbenzamide derivatives. Results indicated that the models have good reliability and predictability. Furthermore, 3D-QSAR contour maps reveal that the bioactivity of antagonists is most affected by hydrogen bond receptors. Molecular docking showed significant hydrogen bonds between key residues (ASN318 and SER103) and OX1 R pocket were determined. Seven new more selective OX1 R antagonists (DS01–DS07) were designed based on the contour maps and molecular docking results. The results of ADME prediction illustrated that the designed compounds had potential druggability and that the selectivity index value was higher than that of template molecule 33 (SI = 265). The molecular dynamics simulation results of compound 33 and DS01 indicated that their stability can be attributed to the hydrogen bonds between the nitrogen atom in the R1 region, the fluorine atom in the R4 region and the protein pockets. We constructed the orexin receptor 2 (OX2 R) QSAR models (CoMFA: q2 = 0.852, r2 = 0.979; CoMSIA: q2 = 0.823, r2 = 0.965) and compared them with models 1. Four potential new skeletal antagonists (Z1–Z4) were screened from the ZINC12 database containing 1,020,679 compounds based on the pharmacophore model. The molecular docking indicated that the screened compounds formed additional hydrogen bond interaction with residues GLN176 and ARG322 of OX1 R. The present study could offer theoretical guidance for future structural optimization, design, and synthesis of more selective OX1 R antagonists. Designed and screened compounds might also be as candidate compounds for related research groups. © 2019 Taiwan Institute of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

1. Introduction Human insomnia disorder is a condition characterized by both nocturnal and diurnal symptoms. It involves a predominant complaint of dissatisfaction with sleep quality or duration and is accompanied by difficulties in initiating sleep at bedtime, frequent or prolonged awakenings, or early-morning awakening with an inability to return to sleep [1]. Statistics show that about 30−35% people worldwide have symptoms of insomnia, making it a threat to human health [2,3]. Traditional insomnia drugs include benzodiazepines, such as temazepam, flurazepam, and triazolam, which can shorten the time of falling asleep, reduce the time, frequency of awakening, and increase the total sleep time. However, use of these drugs can cause dependence, resistance, and memory loss



Corresponding author. E-mail address: [email protected] (Y. Ren).

[4,5]. The later non-benzodiazepines drugs, including zolpidem, eszopiclone, and zaleplon, have certain advantages in terms of tolerance, accompanied with difficulty getting up, sleeping disorders, and residual effect [6,7]. At present, drugs that can fundamentally treat insomnia are lacking. Therefore, new therapeutic drugs for insomnia should be explored. Recent studies have found that human arousal is regulated by a variety of neurotransmitters and neural pathways [8]. OXRs belong to seven-layer membrane G-protein coupled receptor that can emit dense excitatory projections to monoaminergic neurons. Monoaminergic neurons can also emit a back donation on the thalamus and cerebral cortex to maintain an awake state and lead to insomnia [9]. The orexins are divided into two subtypes [10,11], consisting of OX1 R and OX2 R (also known as hypocretin-1 and −2), which are used to develop two types of receptor antagonists. The Biological activity assay was evaluated in CHO (Chinese hamster ovary) cells expressing OX1 R and OX2 R, measuring inhibition of orexin-induced intracellular calcium ([Ca2+ ]i ) release [12]. Hence,

https://doi.org/10.1016/j.jtice.2019.03.018 1876-1070/© 2019 Taiwan Institute of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

2

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

Fig. 1. Representative DORAs, 2-SORAs, and 1-SORAs.

the development of antagonist drugs targeting OXRs has received widespread attention in the field of insomnia drugs [13]. Some representative dual orexin receptor antagonists (DORAs, DORAs are defined as having less than 20-fold selectivity in a cell-based or radioligand binding assay for either OX1 R or OX2 R) [14], selective orexin 2 receptor antagonists (2-SORAs), and selective orexin 1 receptor antagonists (1-SORAs) are shown Fig. 1. Suvorexant is currently the only DORA for the treatment of initial sleep and continuous sleep that was approved by the US FDA in 2014 [15]. Patients who take this medicine may experience drowsiness and sleep hallucinations. It will also have different effects due to differences in sex and weight of the user [16]. Almorexant (DORA) has completed Phase 3 clinical trial in patients with chronic primary insomnia [17]. EMPA, JNJ-10397049, MK-3697, and JNJ-42847922 are 2-SORAs. Some of the 2-SORAs are inefficient, for example, EMPA showed poor uptake in brain in rats and baboon [18]. At present, research on OX1 R antagonists that play an important role in sleepwake effect is still in its infancy. OX1 R antagonists that is really applied to clinical research is currently lacking. Thus, drug development with this target has attracted considerable attention. The first selective OX1 R antagonist that was reported was SB-334867, but its production was discontinued due to lack of oral bioavailability [19]. ACT-335827 was later synthetized as an OX1 R antagonist [20], which addressed the oral problem. The selectivity and pharmacokinetic indicators of ACT-335827 were also defective. A series of pyrazoylethylbenzamide derivatives was designed and synthetize by Futamura (2017) to further explore the pharmacological functions of OX1 R [21]. The pIC50 values of OX1 R and OX2 R are from 6.042–9.115 and 5.836–8.793, respectively, providing a good basis for the models analysis. The structure of the 36 compounds are also diverse. Compound 33 shown in Fig. 2 is a promising selective OX1 R antagonist, which exhibited excellent antagonistic activity against OX1 R, with an IC50 of 2.01 nM and a 265-fold selectivity for OX1 R over OX2 R. However, a systematic theoretical structureactivity relationship study of these compounds remains to be conducted to date. As an effective, rapid, and economical tool, computer-aided drug design has been widely used in the development of various new drugs [22–25]. It not only improves the success rate of

Fig. 2. Structure of compound 33.

new drug development but also provides new strategies for new drug design. Therefore, this study establishes the structure-activity relationship between the structural characteristics and biological activity of pyrazoylethylbenzamide derivatives by building OX1 R QSAR models. Contour maps were used to analyze the biological activity effect of the linked groups in different regions of the antagonist backbone. The binding modes and stability between the OX1 R pocket and the ligand molecules were studied by molecular docking and molecular dynamics simulation. New more selective OX1 R antagonist was designed and predicted for ADME. OX2 R QSAR models were constructed and compared with OX1 R QSAR models. In addition, virtual screening based on the pharmacophore model was performed on the ZINC12 database containing about one million compounds. Flow chart of discovering novel OX1 R antagonists by Computer-Assisted Drug Design is illustrated in Fig. 3. The above systematic research provides important theoretical guidance for the design of more effective 1-SORAs in the future, reduce the blindness of designing such antagonists, ensures the effective synthesis of drugs, and improves the development efficiency. 2. Materials and methods 2.1. Data sets In this study, a total of 36 pyrazoylethylbenzamide derivatives which were designed and synthesized by Futamura were built OX1 R and OX2 R QSAR models. There are two main ways to divide the training set and test set. One of the most widely used

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

3

Fig. 3. Flow chart of discovering novel OX1 R antagonists by computer-assisted drug design.

method is a mere random selection. Another frequently used approach is based on the activity sampling. The whole range of activities is divided into bins, and compounds belonging to each bin are randomly assigned to the training or test sets [26]. Type 2 is more suitable for our study and selected in our work. Training set composed of 29 compounds for models construction; the other is a test set consisting of 7 compounds. The purpose of the test set is to verify the reliability of the training set model. The IC50 value of the nanomolar bioactivity of the compound is known, and the IC50 value is converted to the corresponding pIC50 value (pIC50 = −logIC50 ) for analysis of the model dependent variable. The molecular structure, experimental, and predicted CoMFA, CoMSIA, and HQSAR pIC50 values of the OX1 R models and OX2 R models were indicated in Tables S1 and Table S2, respectively. 2.2. Molecular construction and alignment The three-dimensional structures of the 36 pyrazoylethylbenzamide derivatives were constructed in the Sketch module of Sybyl-X 2.0 software package (Tripos Inc. St. Louis, USA) and the

original conformations of all molecules were decided according to the conformation of the Suvorexant extracted from the proteins (OX1 R PDB cod: 4ZJ8; OX2 R PDB cod: 4S0V) [27,28]. Then, each of the constructed molecule was submitted to energy optimization. The purpose of energy optimization is to obtain a stable molecular bioactive conformation that binds to the receptor. The Powell method was optimized here and the Tripos force field of the Gasteiger–Hückel atomic charge was used. The maximum number of iterations was 10 0 0 and the energy convergence gradient was set to 0.005 kcal/(mol A). The other parameters were set to the system default values. The alignment operation is not only a fundamental step to establish reliable 3D-QSAR models but also the quality of the alignment can directly affect the predictive ability of the models. Compound 12 has the highest activity (IC50 = 0.767 nM) for OX1 R was selected as a template molecule for alignment. Fig. 4(A) is the chemical structure of the template molecule 12 and purple atoms represent common skeleton. Fig. 4(B) illustrates the alignment result based on the common skeleton. Compound 11 which has the best activity (IC50 = 1.61 nM) for OX2 R was chose as a

4

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

In addition to the aforementioned internal validation, the reliability of the models can be further verified by external validation. External validation predicted correlation coefficient (r2 pred ) has the same validation ability as the statistical parameters calculated by the aforementioned Partial Least Square (PLS) analysis. The predicted correlation coefficient (r2 pred ) that higher than 0.5 may be considered as an indicator of effective external predictability [34]. r2 pred equation: Fig. 4. (A) The chemical structure of template molecule 12 and common skeleton marked in fuchsia. (B) The alignment results of training set.

template molecule for alignment. Fig. S1A is the chemical structure of the template molecule 11 and the purple atoms represent common skeleton. Alignment result based on the common skeleton is shown in Fig. S1B. 2.3. CoMFA, CoMSIA, and HQSAR CoMFA and CoMSIA are two widely used three-dimensional quantitative structure-activity relationship tools [29,30]. Two models place compounds in a 2 A˚ spatial lattice, which links the three-dimensional structure, biological activity data, and the potential energy of the molecule, providing valuable information for the molecular modification. The CoMFA model is based on the Lennard–Jones potential energy function and Coulomb’s law, which simulates sp3 hybridized carbon atoms and +1 valent hydrogen atoms as probe atoms, respectively. Due to some problems of the algorithm itself, the cut-off value is specifically defined as 30 kcal/mol. However, CoMSIA is completely different from CoMFA, which adopts the distance-associated Gaussian function. The CoMSIA model is calculated using a sp1 hybridized carbon atoms as probe atoms. In addition, the CoMSIA models no longer limits the steric field and the electrostatic field but also covers the hydrophobic field, hydrogen bond donor field, and hydrogen bond receptor field, which also make up some shortcomings of CoMFA itself. In a short, both CoMFA and CoMSIA provide a better explanation for the three-dimensional structure activity relationship of compounds. HQSAR is a fragment-based method that need not for compounds 3D structure and molecular alignment [31]. Evaluation studies have shown that HQSAR has better predictive power than those complex 3D-QSAR techniques. HQSAR model can be affected by a number of parameters concerning hologram generation: hologram length, fragment size, and fragment distinction. The information in each fragment is defined by fragment distinction parameters: atoms (A), bonds (b), connections (C), hydrogen atoms (H), chirality (Ch), and donor acceptor (DA) [32]. Each fragment is assigned a specific positive integer in the range of 0–231 via a cyclic redundancy check algorithm.

r 2pred = 1 −

P RESS SD

(1)

In the above formula, PRESS means the sum of squares of the deviation between the experimental and predictive pIC50 values of the test set molecule, and SD is the sum of squares of the deviation between the experimental pIC50 value of test set molecule and the average experimental pIC50 value of all training set molecules. According to Roy and Tropsha [35,36], the other four statistical parameters can more strongly verify the reliability of the model, and these parameters need to satisfy the following conditions:

R2 > 0.6, 0.85 ≤ K ≤ 1.15,







R2 − R20 /R2 < 0.1, R2m > 0.5.

R2 is the squared correlation coefficient between the predicted and experimental values of the activity. Equation is as follows:

    2     Y t est − Y t est Y pred ( test ) − Y pred ( test )  R2 =          Y t est − Y t est 2  Y pred (test ) − Y pred (test ) 2  (2) K is the slopes of regression lines through the origin for fitting to experimental and predicted data. Equation is as follows:

 K=

(Ytest ∗ Ypred (test ) )  (Ypred (test ) )2

(3)

R20 is the squared correlation coefficient obtained by using predicted versus observed activities. Equation is as follows:

 R20 = 1 − 

(Ypred (test ) − kYpred (test ) )2 2

(Ypred (test ) − Y pred (test ) )

(4)

R2 m is another criteria proposed by Roy to investigate the external predictability of QSAR models. Equation is as follows:

R2m = R2 ∗ (1 −



R2 − R20 )

(5)

Where Ytest , Y test, Ypred(test), and Y pred (test ) are experimental, average, predicted and the predicted average activities in the test set, respectively.

2.4. Model validation 2.5. Molecular docking Model validation is a key step in building any QSAR model. Partial least square (PLS) analysis is a multivariate regression analysis method based on linear transformation ideas, which combines the potential fields and bioactivity data involved in CoMFA and CoMSIA [33]. Leave-one-out (LOO) cross-validation method was utilized to obtain the cross-validation correlation coefficient (q2 ) value and optimum number of components (ONC). Then, on the basis of optimum number of components, non-cross-validation was performed, and the non-cross-validation correlation coefficient (r2 ) value, standard error estimate (SEE), and F-test value were obtained. In addition, a p-value of 0.05 was used to determine the Fisher test’s significance in the model building process by SPSS (version 22.0). These significant statistical parameters reflect the reliability of the built models.

Molecular docking is an effective strategy to study the microscopic interaction between ligands and receptors, which can give us significant information about them [37]. Molecular docking simulations were operated on the Surflex-Dock module of the Sybyl-X 2.0 software [38]. Two crystal structures of OX1 R (PDB code: 4ZJ8) and OX2 R (PDB code: 4S0V) can be obtained from Protein Data Bank (https://www.rcsb.org/) [39]. There are two main steps in docking simulations. The first step is the preparation of two OXRs. Firstly, all water molecules of OXRs were removed from the original 4ZJ8 and 4S0V. Side chains and termini chains were repaired after analyzing the protein. The hydrogen atoms and the essential charges were then added to the protein. Secondly, compounds 33, 36, and the newly design compound DS01 (shown in Table 3) were

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17 Table 1 The statistical values of pharmacophore models of GALAHAD. NO.

Specificity

N_Hits

Feats

Pareto

Model-001 Model-002 Model-003 Model-004 Model-005

5.620 4.970 2.715 2.543 2.540

20 18 18 5 10

8 10 10 11 11

0 0 0 0 0

docked into OX1 R and OX2 R, respectively. Finally, ten different conformations for each compound were obtained. 2.6. Molecular dynamics (MD) simulation Molecular dynamics simulation were conducted in the Amber 14.0 software package [40]. Compound 33 and DS01 were utilized to perform molecular dynamics simulation. Before the simulation, molecular 33 and molecular DS01 were merged into prepared OX1 R pocket, respectively. In the tleap module, the OXRs and ligand were performed for FF99SB force field and a generally general amber force field (gaff), respectively. The ligand-protein complexes were simulated in a hexahedral TIP3P water box of size 78 A˚ ∗ 78 A ˚ ∗ 78 A˚ . The system charge was then neutralized by the addition of Cl− . The minimization, heating, balance and production procedures were then performed on the dynamics simulation environment. Finally, a 10-ns molecular dynamics simulation was conducted in the Amber 14.0 software package. 2.7. ADME prediction Pharmacokinetics reflects the basic rules of absorption, distribution, metabolism, and excretion of drugs in the human body. According to statistics, most of the candidate drugs failed to pass the clinical trial due to poor pharmacokinetics. Hence, predicting ADME properties of drugs has become an indispensable part of developing new drugs [41]. In order to evaluate the pharmacokinetics of our newly designed selective OX1 R antagonists, the SwissADME web tool was adopted through the website http://www. swissadme.ch [42]. The parameters included in the analysis are logP (octanol-water partition coefficient), TPAS (Topological polar surface area), log S (water solubility), the number of rotatable bonds, HIA (Human gastrointestinal absorption), BBB (blood-brain barrier), CYP2C9 (Cytochrome P450-2C9), and synthetic accessibility. Predicting ADME properties of compounds will provide theoretical supporting. 2.8. Virtual screening In the process of discovering new drugs, the application of virtual screening not only promptly accelerate the discovery of lead compounds, but also save the huge cost of blind experimental screening [43]. Virtual screening based on ligands is a wellrecognized solution and has been applied to experiments [44]. Two preparations were required of virtual screening work, including pharmacophore construction and database preparation. In this study, five pharmacophore models were generated by five different groups of compounds (Table S3) with the genetic algorithm of database (GALAHAD) module in the SYBYL-2.0 software. Statistical values of five pharmacophore models were shown in Table 1. Each model is given a Pareto rank 0, which indicates that no one model is superior to any other. The best pharmacophore model (Model-001,Specificity = 5.620) was selected to screen the ZINC12 database [45], containing a total of 1020,679 molecules by using UNITY Flex search program.

5

ZINC12 database can obtained by accessing online webpage (http: //zinc.docking.org/). Compounds with more than 40 QFIT values were further studied by docking them into OX1 R by Surflex-Dock. 3. Results and discussions 3.1. Discussion of CoMFA, CoMSIA, and HQSAR results The q2 of excellent 3D-QSAR models should be higher than 0.5, whereas r2 should be greater than 0.9, along with low SEE and great F values. 3D-QSAR models parameters of OX1 R were summarized in Table 2. There are many different models, including ED, EHD, SEHDA, DA, and SD of CoMSIA models. S, E, H, D, and A represent steric field, electrostatic field, hydrophobic field, hydrogen bond donor field, and hydrogen bond acceptor field, respectively. For example, CoMSIA-ED indicates the effect of electrostatic field and hydrogen bond donor field on the biological activity of the antagonists. The same is true in other models. The q2 of all the models were over 0.7 except CoMSIA-ED (q2 = 0.494) and CoMSIA-EHD (q2 = 0.692), suggesting that nearly all of the established models were reliable. The corresponding field contributions were 67.1% of S and 32.9% of E. Meanwhile, it’s common to use combinations of S, E, H, D, and A to generate the CoMSIA model. We found that the electrostatic field contribution rate (24.1%) of is much larger than the electrostatic field contribution rate (7.3%) of CoM SIA-SEHDA model. As the number of fields increases, the contribution rate of the steric field has a tendency to decrease in all CoMSIA models. Among ED, DH, and DA of CoMSIA models, the contribution rate of D was weaker than that of other field. Electrostatic field, hydrophobic field, and hydrogen bond acceptor field made the most influence on the bioactivity. Finally, CoMSIA-SEHDA model was chosen as the best 3D-QSAR models. In the OX2 R models, Table S4 showed that the q2 of all the models were more than 0.7 except CoMSIA-SE (q2 = 0.690), CoMSIA-SA (q2 = 0.578), and CoMSIA-DA (q2 = 0.639), indicating that nearly all of the established models were reliable. The effects of steric field and the electrostatic field were basically same in the CoMFA model. In the CoMSIA, the electrostatic field contribution was significantly larger than the steric field. In the same way, CoMSIA-SEHDA was selected as the best CoMSIA model at last. To test the predictive power of OX1 R and OX2 R models, 7 test sets were selected to assess the predictive abilities. The predictive correlations (r2 pred ) of CoMFA-SE were 0.787 and 0.846 and CoMSIA-SEHDA were 0.817 and 0.797, indicating the good external predictive capacity. Scatter plots of experimental and predicted pIC50 values of training sets (marked in black) and test set (marked in red) for the two 3D-QSAR models were presented in Figs. 5 and S2. The distribution of all points has a significant linear regression law, indicating that there were small differences between experimental activity and predictive activity. In a summary, OX1 R and OX2 R 3D-QSAR models we established had good predictability and reliability. In addition, p value of 0.05-quantile of Fisher’s distribution were summarized in Table S5. The p values of OX1 R and OX2 R QSAR models were all over 0.05, which reveals that there were no significance level between experimental and predicted values. Tables 2 and S4 summarized the two optimal HQSAR model statistical parameters of OX1 R and OX2 R. The best HQSAR model of OX1 R was generated using atoms, Bonds, Connections, Hydrogen Atoms, Chirality, and Donor Acceptor as fragment distinction and 4–7 as the fragment size at the hologram length (HL) of 83. The HQSAR model uses a color code to discriminate the main atomic contributions to the activity. The molecule is color coded to reflect the individual atomic contributions. The colors at the red end of the spectrum (red, red orange, and orange) reflect poor (or negative) contributions, while colors at the green end (yellow,

6

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17 Table 2 CoMFA, CoMSIA, and HQSAR statistical analysis results of the OX1 R. Parameters

CoMFA

CoMSIA

HQSAR

q2 r2 SEE F ONC r2 pred R2 k R2 0 (R2 − R0 2 )/R2 R2 m

0.788 0.988 0.090 255.847 7 0.787 0.888 1.0 0 0 1.0 0 0 −0.126 0.590

0.783 0.953 0.185 50.849 8 0.817 0.934 0.998 0.998 0.069 0.696

0.512 0.927 – – 6 0.917 0.935 0.998 0.981 −0.049 0.734

PLS statistics

q2

r2

SEE

F

ONC

Contribution (%) S

E

H

D

A

COMFA-SE COMSIA-SE COMSIA-SH COMSIA-SD COMSIA-SA COMSIA-EH COMSIA-ED COMSIA-EA COMSIA-HD COMSIA-HA COMSIA-DA COMSIA-SEH COMSIA-SED COMSIA-SEA COMSIA-SHD COMSIA-SHA COMSIA-SDA COMSIA-EHD COMSIA-EHA COMSIA-EDA COMSIA-HDA COMSIA-SEHD COMSIA-SEHA COMSIA-SEDA COMSIA-SHDA COMSIA-EHDA COMSIA-SEHDA

0.788 0.822 0.802 0.754 0.766 0.736 0.494 0.705 0.744 0.804 0.769 0.785 0.806 0.782 0.802 0.809 0.802 0.692 0.770 0.784 0.788 0.771 0.790 0.789 0.805 0.779 0.783

0.988 0.977 0.960 0.941 0.969 0.971 0.944 0.962 0.963 0.967 0.931 0.973 0.975 0.970 0.967 0.969 0.974 0.967 0.968 0.935 0.948 0.973 0.969 0.944 0.973 0.949 0.953

0.090 0.136 0.172 0.203 0.158 0.154 0.212 0.176 0.169 0.165 0.219 0.149 0.143 0.156 0.158 0.158 0.146 0.164 0.162 0.208 0.195 0.149 0.158 0.197 0.147 0.193 0.185

255.847 77.654 59.388 47.589 56.835 59.618 30.556 45.568 54.464 52.127 40.302 63.949 69.541 58.201 62.764 56.994 66.582 52.389 53.810 52.322 45.690 63.974 56.650 50.549 65.620 46.297 50.849

7 10 8 7 10 10 10 10 9 10 7 10 10 10 9 10 10 10 10 6 8 10 10 7 10 8 8

67.1 31.4 25.0 38.7 36.4 – – – – – – 18.4 25.6 22.8 19.3 16.6 21.4 – – – – 15.5 13.9 9.1 11.4 – 7.3

32.9 68.6 – – – 45.8 62.9 58.5 – – – 40.0 48.7 42.0 – – – 36.7 31.8 32.8 – 32.5 28.1 30.9 – 25.6 24.1

– – 75.0 – – 54.2 – – 71.2 61.4 – 41.6 – – 57.7 46.8 – 44.9 38.9 – 41.5 36.7 30.3 – 37.9 22.4 21.1

– – – 61.3 – – 37.1 – 28.8 – 27.0 – 25.7 – 23.0 – 26.7 18.4 – 23.2 17.3 15.3 – 20.5 14.1 18.5 16.9

– – – – 63.6 – – 41.5 – 38.6 73.0 – – 35.2 – 36.6 51.9 – 29.3 44.1 41.2 – 27.7 39.4 36.7 33.5 30.8

Fig. 5. Plots of experimental and predicted pIC50 values for the total set in the CoMFA (A), CoMSIA (B), and HQSAR (C) models of the OX1 R.

green blue, and green) reflect favorable (positive) contributions. Atoms with intermediate contributions are colored white. The HQSAR contribution map of compound 12 and 11 were illustrated in Figs. 6 and S3. 3.2. 3D-QSAR contour maps analysis of OX1 R and OX2 R models Analyzing the CoMFA and CoMSIA models, some physicochemical properties of the potential field around the compound itself

can be quantitatively reflected by the contour maps. Here, the relationship of five potential fields under the three-dimensional quantitative structure-activity relationship analysis was discussed. In the steric field, a green color block indicates that introduction of a bulky group in this region increases the biological activity of the antagonists, whereas yellow contour maps mean introducing a large group can reduce potency. In the electrostatic field, the red contour maps with electronegative groups were favorable to the activity, whereas the blue contour maps correspond to the

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

Fig. 6. HQSAR contribution map of compound 12. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

regions where electropositive groups were favorite. In the hydrophobic field, the yellow contour maps mean that the hydrophobic groups introduced in this area are useful to improve the activity of the antagonists, whereas the silver contour maps indicate that the hydrophilic groups introduced in this area are useful to improve the activity of the compounds. In the hydrophobic field, yellow indicates that the introduction of a hydrophobic group is advantageous to increase the potency of the antagonists, and silver indicates that the introduction of a hydrophilic group is more advantageous. In the hydrogen bond field, blue–green indicates that the introduction of a hydrogen bond donor groups is helpful to increase the activity of the antagonists, while purple is opposite. The magenta color indicates that introduction of a hydrogen bond acceptor groups is advantageous to enhance the activity of the compounds, while red is the opposite. Fig. 7(A) and (E) is OX1 R CoMFA and CoMSIA steric field contour maps, respectively. In Fig. 7(A), it can be obviously seen that many yellow contour maps were displayed in the vicinity of the R1 and R4 regions, whereas a relatively large green contour map existed at the linkage moiety. Through the CoMFA steric field color maps, the size of the compounds bulk tend to be relatively small. Moreover, in Fig. 7(E), the molecular 33 was surrounded by a huge yellow contour map and a relatively small green contour map in the linkage moiety, which is consistent with the CoMFA model analysis. Adding some relatively small groups in the R1 region will increase the activity. The bulk of the triazole group is smaller than the pyrimidine group so the pIC50 values of the compounds will be ranked as follow: 12 > 19, 20 > 21, 23 > 24, 26 > 27, 18 > 28, 29 > 30, 35 > 36. At the linkage moiety, adding some relatively large groups will enhance the activity of the molecule. The group bulk sequence decreasingly is: ethyl, methyl, hydrogen atom. For example: 05 > 03 > 01, 06 > 04 > 02, 11 > 09 > 07, 12 > 10 > 08, and 07–12 are greater than 01–06, respectively. However, there is an exception, such as 09 > 14 > 15. The R3 groups of three compounds were linked to a methyl group, an ethyl group, and an isopropyl group, respectively. It is demonstrated that adding a large group at the linkage moiety increases the biological activity but it is improper when the group is relatively large. It may be due to the fact that the size of the active site pocket is relatively fixed and it is not suitable for a relatively big compound, resulting in decreased activity. Fig. 7(B) and (F) are OX2 R CoMFA and CoMSIA steric field contour maps, respectively. In Fig. 7(B), it can be clearly found that the contour maps distribution of R1 region and linkage moiety are basically the same as those in Fig. 7(A). It is worth noting that the yellow contour maps in the R4 region are significantly reduced compared to Fig. 7(A), indicating that it is not appropriate to introduce small groups. Fig. 7(B) is the CoMFA electrostatic field contour maps of the OX1 R models. At the R3 region, there was a medium-sized red contour map, indicating that the electronegative group is favorable. The electronegativity of the methyl is greater than the hydrogen atom so the biological activity of the compounds 07–12

7

are greater than 01–06, respectively, which is consistent with steric field analysis of the previous OX1 R models. A larger blue contour map was above the R4 region, indicating that having electronegative substituents here is disadvantageous for increasing the activity. The electronegativity of the nitrogen atom is greater than the carbon atom, making the compounds attached to the pyridine ring inferior to the compounds attached to the benzene ring. For example, 20 > 17, 21 > 19, 23 > 22, 29 > 18, 30 > 28, 32 > 31. Fig. 7(G) is contour maps distribution of CoMSIA electrostatic field, which is basically consistent with the CoMFA model. Fig. 7(D) is the CoMFA electrostatic field contour maps of the OX2 R models. It can be found that there were a red contour map and a blue contour map above the R3 and R4 region, respectively, and the bioactivity sequence of the previous compounds is consistent with the model. It is worth noting that there were two medium-sized red contour maps on the inside of the R4 region, which were not found in Model 1, indicating that it is advantageous to add an electronegative group here. Fig. 7(H) is contour maps distribution of CoMSIA electrostatic field, which is basically consistent with the CoMFA model. Fig. 8(A) shows the CoMSIA hydrophobic field contour maps of the OX1 R models. A huge yellow contour map was embedded in the R4 region, indicating that the hydrophobic group is favorable. The hydrophobicity of the benzene ring is greater than the pyridine ring and the degree of extension of the benzene ring to the yellow contour map is correspondingly greater than the pyridine ring. So the compound activity of 20, 21, 29, 30 are more potent than the 17, 19, 18, 28, respectively. However, there are two exceptions, such as (22–23) and (31–32). The benzene rings of these two groups do not extend to the yellow contour maps as much as the pyridine ring. Interestingly, the R1 region of the two group compounds are attached to fluorine atom at the benzene ring 5-position, which is the biggest difference from the previous four groups (the methyl group attached at the 5-position). This is consistent with the previous analysis of electrostatic field activity. There is a silver contour map above the tetrazole ring in the R4 region, indicating that the hydrophilic group is favorable here. Fig. 8(B) illustrate the CoMSIA hydrophobic field contour maps of the OX2 R models. It can be found that an enormous yellow contour map was embedded in the R1 region and the R4 region. Compared with the distribution of the hydrophobic field contour maps in the OX1 R models, there was a yellow contour map in the R1 region and a silver contour map was missing in the tetrazole ring. It is an effective strategy to consider modifying the two different regions for new selective OX1 R compounds. Fig. 8(C) is the CoMSIA hydrogen bond field contour maps of the OX1 R models. A medium-sized magenta contour map was located above the pyrimidine ring 2-position in the R1 region, indicating that the hydrogen bond acceptor is favorable. This is a possible reason why these compounds which contain a triazole ring (compounds: 23, 26, 28, 29, 35) were more potent than those that connect pyrimidine ring(compounds: 19, 21, 24, 27, 28, 30, 36), which is consistent with the conclusions obtained from the previous steric field. In addition, a relatively large red contour map was above the fluorobenzene ring in the R4 region, illustrating that the hydrogen bond acceptor is disadvantageous here. There was a large red contour map above the benzene ring 2-position and the activity of the pyridine ring compounds (17, 19, 18, 28, 22, 31) is lower than those containing benzene ring compounds (20, 21, 29, 30, 23, 32), respectively. It is consistent with previous electrostatic field and hydrophobic field analysis results. Fig. 8(D) is the CoMSIA hydrogen bond contour maps distribution of the OX2 R models. There was a medium-sized magenta contour map under the pyrimidine ring 6-position and a medium-sized red above the fluorobenzene 6-positionn in the R1 region which were not observed in the model 1.

8

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

Fig. 7. OX1 R CoMFA and CoMSIA steric fields (A and E) and electrostatic field (C and G) contour maps. OX2 R CoMFA and CoMSIA steric fields (B and F) and electrostatic fields (D and H) contour maps. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 8. OX1 R CoMSIA hydrophobic field (A) and hydrogen bond field (C) contour maps. OX2 R CoMSIA hydrophobic field (B) and hydrogen bond field (D) contour maps.

3.3. Design of novel OX1 R compounds By analyzing the three-dimensional quantitative structureactivity relationship between the OX1 R and OX2 R, novel selective OX1 R antagonists via compound 33 as a template molecule were designed. Structure-activity relationship (SAR) information of OX1 R (A) and OX2 R (B) are illustrated in Fig. 9. By comparing and analyzing the distribution of contour maps between OX1 R and OX2 R, the largest differences between them were R1 region and R4 region. In the electrostatic field of the OX2 R CoMFA model, a large blue contour map extended to the fluorobenzene ring 6-position in the R1 region (there was no such a contour map in the same region of the OX1 R models), indicating electropositive substituents are favorite. Hence, adding an elec-

tronegative group fluorine atom here can reduce the biological activity of the antagonist to the OX2 R. In the CoMSIA hydrophobic field of the OX2 R models, there was a favorable yellow hydrophobic contour. So adding a hydrophilic group carboxyl can reduce the antagonist to the OX2 R bioactivity. In the CoMSIA hydrogen bond field of the OX2 R models, there was a medium-sized magenta contour map at the 6-position nitrogen atom of the pyrimidine ring at the R1 region (the OX1 R models did not have this pyrimidine in the same region), indicating hydrogen bond acceptors are favorite. Therefore, modifying the pyrimidine ring to pyridine would reduce the bioactivity of the antagonists to the OX2 R. In the CoMSIA hydrophobic field and hydrogen bond field of the OX1 R models, there was a small silver contour maps near the tetrazole ring in the R4 region and a medium-sized red contour maps at the

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

9

Fig. 9. Structure-activity relationship (SAR) information of OX1 R (A) and OX2 R (B).

Fig. 10. Re-docking results of the crystal structure (B, colored in green) into the pocket of OX1 R (A) and the cognate ligand (B colored in cyan) and the crystal structure (D, colored in orange) into the pocket of OX2 R 2 (C) and the cognate ligand (D colored in cyan). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

fluorobenzene 6-position. So adding Hydrophilic hydrogen bonds donor group hydroxyl and carboxyl groups would enhance the selectivity for OX1 R over OX2 R. Here, compounds DS01–DS07 illustrated in Table 3 were designed via compound 33 as a template molecule. It found that all the designed compounds were more selective than the compound 33 except for the compound DS06. The strategy we adopted was to reduce the bioactivity of the compounds to the OX2 R. 3.4. Molecular docking analysis Molecular docking is a common mean of simulating proteinligand interactions by computer software. By docking, important information about one or more compounds at a particular protein site can be provided, including the overall geometric conformational orientation of the compounds and microscopic interactions within the cavity. In order to illustrate the accuracy of OX1 R and OX2 R pockets, the target ligand (Suvorexant) was re-docked into the crystal structures of OX1 R protein (PDB code: 4ZJ8) and OX2 R protein (PDB code: 4S0V). The re-docking results were illustrated

in Fig. 10. It was found that the original ligand and re-docked ligand maintain a high degree of uniformity in the overall conformational orientation. In this study, the corresponding similarity of the re-docking results were 0.904 and 0.872, respectively, which indicated that the docking method was rational and adoptable. The non-bonding interactions between the molecule and the active site residues stabilize the compound in the cavity. Hence, the two different active site pockets described above can be further used for the next docking analysis. Among the pyrazoylethylbenzamide derivatives, the compound 33 exhibited 265-fold selectivity for OX1 R over OX2 R. In order to gain a deeper understanding of the microscopic interaction patterns of compound 33 in the two active sites, molecular 33 was docked into two subtype pockets to analyze differences. Table 4 shows significant amino acid residues which present in ˚ The main differences bethe OX1 R and OX2 R pockets within 3 A. tween the two pockets were a substitution of THR135 for ALA127 and THR111 for SER103 in OX2 R. VAL347 and VAL106 were distributed obliquely below the pyridine ring of Suvorexant R1 region and obliquely above the SER103, respectively. However, TYR354

10

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17 Table 3 Structures of designed molecules and their predicted pIC50 based on the OX1 R and OX2 R models.

NO.

R1

R4

Predicted pIC50 (4ZJ8)

Predicted pIC50 (4S0V)

CoMFA

CoMSIA

CoMFA

CoMSIA

CoMFA Residual

CoMSIA residual

33

8.671

8.638

6.242

6.086

2.429

2.552

DS01

8.274

8.549

4.969

5.558

3.305

2.991

DS02

7.916

7.705

5.352

5.166

2.564

2.639

DS03

7.886

8.116

5.398

5.504

2.488

2.612

DS04

8.041

7.854

5.414

5.215

2.627

2.639

DS05

8.022

8.139

5.281

5.05

2.741

3.089

DS06

7.916

7.469

5.244

5.013

2.672

2.456

DS07

7.95

7.705

5.293

4.846

2.657

2.859

and MET191 were distributed directly below the THR111 of OX2 R pocket and the tetrazole ring of Suvorexant R4 region, respectively. The aforementioned four groups different amino acid residues were one of the biggest differences of OX1 R over OX2 R pockets. Fig. 11(A) and (C) was the pocket maps of compound 33 docked into OX1 R and OX2 R. Fig. 11(B) and (D) was corresponding 2D nonbonded interaction maps. Compound 33, like Suvorexant, accepts a horseshoe-like conformation in both pockets [46]. In the OX2 R pocket, there was an F-shaped π –π stacking between fluorobenzene ring and tetrazole ring. Fluorobenzene in the R1 region not only forms a π -π stacking interaction with the tetrazole in the R4 region, but also has another π -π stacking interaction with the fluorobenzene in this region, which means the interaction with OX1 R pockets was stronger. It was found that the R4 region in the two

pockets had a certain degree decline relatives to original ligand. The key residue THR111 of OX2 R and SER103 of OX1 R interacted with the antagonist 33 by hydrogen bonding, respectively. The hydrogen bond distances observed were 2.889 A˚ (THR111-NH • • • F) and 2.811 A˚ (SER103-NH • • • F), which indicated that molecule 33 combined with OX1 R pocket was tighter than the OX2 R pocket. Moreover, the new chain skeleton linked R1 and R4 regions instead of a ship skeleton like Suvorexant, which means that the interaction between the molecule 33 and R4 region was relatively weak than Suvorexant. In addition, residue ASN318 of OX1 R and residue ASN324 of OX2 R formed significant hydrogen bond with molecular 33. The Jie Yin group has explained the importance of the hydrogen bond here. It is worth noting that the pyrimidine ring of R1 region was rotated to different degrees angle in both pockets. In

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

11

Fig. 11. Docking results and 2D non-bonded interaction maps of selected compounds 33 in the binding site of OX1 R pocket (A, B) and OX2 R pocket (C, D).

Table 4 Two subtypes of protein ligand molecules Important ˚ amino acid residues within 3 A. PDB

4ZJ8

4S0V

Released Method Resolution 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

2016–03–09 X-ray diffraction 2.751A˚ HIS344 ASN318 ILE314 TYR311 HIS216 PHE219 VAL130 GLN179 PRO123 GLN126 ILE122 CYS99 TRP112 ALA102 ALA127 SER103 VAL347 VAL106

2015–01–14 X-ray diffraction 2.5 A˚ HIS350 ASN324 ILE320 TYR317 HIS224 PHE227 VAL138 GLN187 PRO131 GLN134 ILE130 CYS107 TRP120 ALA110 THR135 THR111 TYR354 MET191

the OX1 R pocket, the pyrimidine ring formed a 2.163 A˚ hydrogen bond with the residue ASN318 (ASN318-NH • • • N) which stabilized molecular 33 with pocket. The pyrimidine ring rotated outward 45° geometrical angles. However, the pyrimidine ring rotated outward almost 90° geometrical angles due to lacking a hydrogen bond with key residue ASN324 in OX2 R pocket. This different spatial conformation exhibited by the pyrimidine ring of compound 33 in different pockets was likely to be one of the main reasons for 265-fold selectivity for OX1 R over OX2 R. Moreover, in the pre-

vious CoMSIA model analysis, the contribution of hydrogen bond acceptor accounted for 20.8% and 22.9% in two models, respectively, indicating that the hydrogen bond receptor has an advantage in the OX1 R models. In addition, the hydrophobic interaction between the compound 33 and significant residues of OX1 R such as ALA102, PRO123, ALA127, HIS216, PHE219, ILE314 and VAL347 further enhanced the antagonist bioactivity. So molecule 33 was more stable to OX1 R than to OX2 R. Compound 36 had 213-fold selectivity for OX1 R over OX2 R by replacement of the 5-methylpyridine with 5-fluorobenzene, which exhibited better antagonistic activity against OX1 R with an IC50 value of 2.01 nM. Comparing compound 33, IC50 of molecular 36 increased by 4.14 nM and 778 nM for OX1 R and OX2 R, respectively. Docking results of molecular 36 into OX1 R and OX2 R pockets were shown in Fig. 12. In the OX1 R pocket, compound 36 still adopted horseshoe-shaped bioactive conformation. Comparing compound 33, the fluorobenzene of the R4 and pyrimidine ring of the R1 region of molecule 36 did not form hydrogen bonds with residue SER103 and ASN318, respectively. It was noteworthy that there ˚ was a stronger 2.096 A-hydrogen bond interaction between residue ASN318 and amide (ASN318-NH • • • O=) while hydrogen bond ˚ Although the length of compound 36 with ASN318 was 2.496 A. number of hydrogen bonds reduced, the stability of the hydrogen bond formed by the molecule 36 and the residue ASN318 was improved. In addition, residue HIS344, pyridine of R1 region, tetrazole ring and benzene ring of R4 region and residue TRP112 formed two F-shaped and one T-shaped π –π stacking, which illustrated that the multi-conjugate effects formed by the groups further consolidated the degree of binding of the ligand to the receptor. However, molecule 36 did not form additional two hydrogen bonds with OX1 R, which was likely to be one of the main reasons why the latter activity value was lower than the former. The binding conformation and interaction of the molecule 36 in the OX2 R were

12

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

Fig. 12. Docking results and 2D non-bonded interaction maps of selected compounds 36 in the binding site of OX1 R pocket (A, B) and OX2 R pocket (C, D).

shown in Fig. 12(C) and (D). The biggest difference was that molecular 36 has deviated from the horseshoe-shaped bioactive conformation, which may be one of the significant factors for the sudden decline in the activity of the compound 36. From Jie Yin’s viewpoint, the main differences between the two pockets are a substitution of THR111 for ALA127 and THR135 for SER103 in OX2 R, both of which result in a 5% decrease (32 A˚ 3 ) in the volume of the pocket relative to OX1 R. The molecule 36 must bind to the pocket 2 to resolve the clash, which lets compound 40 undergo dramatic changes in conformation and have poor activity. Among 7 designed compounds, the predicted activity of compound DS01 was similar to compound 33 by the OX1 R CoMFA and CoMSIA models. While predicted activity was much lower than compound 33 by the OX2 R CoMFA and CoMSIA models, which improving the selectivity for OX1 R over OX2 R. Docking results and 2D non-bonded interaction maps of designed compounds 01 in the binding site of the OX1 R pocket (A, B) and OX2 R pocket (C, D) were shown in Fig. 13. Hydrogen bond between the compound DS01 and the key residue ASN318, SER103 of orexin pocket 1 was still existing. It was noteworthy that there were three stable hydrogen bonds between the amide oxygen atom, pyrimidine ring nitrogen atom and the residue ASN318, which enhanced the stability of molecule DS01. Comparing to compound 33, the hydrogen bond distance between the amide oxygen atom and residue ASN318 was 1.933 A˚ in OX2 R pocket, which illustrated the interaction was weaker.

3.5. Molecular dynamics (MD) simulation analysis Molecular dynamics simulation can effectively describe the trajectory of each atom in the molecular system over a period of time. To further explore the true binding patterns between OX1 R pocket and the antagonists and verify the accuracy of the docking results,

two complexes (4ZJ8-33 and 4ZJ8-DS01) were performed a 10-ns molecular dynamics simulation by Amber 14.0 software package. The-root-mean-square deviation (RMSD) of the structural skeleton is a data that can describe the differences between the initial structure and the final structure in the protein versus time, which has been an important value for judging whether the system reached balance. The RMSD trajectory of the two complexes (4ZJ8-33, blue; 4ZJ8-DS01, red) were shown in Fig. 14. It can be seen that the RMSD values of the two complexes were balanced after a 10-ns molecular dynamics simulation, which gave about 4 A˚ of 4ZJ8-33 and 3.5 A˚ of 4ZJ8-DS01, respectively. It is worth noting that the RMSD values of the two complexes systems has been in a steady rising phase during the first 3-ns simulation process. At ˚ 3 ns, the RMSD values of the two were basically ame, about 3.5 A. However, after 3 ns, the RMSD value of the 4ZJ8-33 complex system gradually increased and finally stabilized at around 4 A˚ while ˚ 4ZJ8-DS01 complex system was stable at around 3.5 A. The initial starting structures with poses from the final frame of complexes (4ZJ8-33 and 4ZJ8-DS01) were shown in Fig. 15. In the OX1 R pocket, the starting and ending ligand DS01 positions align well (Fig. 15(B)) while the pyrimidine ring of molecule 33 rose. Comparing molecule 33, compound DS01 is more similar to the initial conformation and its RMSD value is relatively more stable. It can also be seen that the conformation variation of the two ligands in OX1 R pocket changed very little by the dynamic alignment of the two complexes, which indirectly illustrated the reliability of the molecular docking that we simulated earlier. By analyzing dynamics simulation results, we can find that some important residues, such as ASN318, SER103, ILE314, VAL347, still play important role in the interaction between antagonist DS01 and OX1 R protein. As with the previous molecular docking results, the hydrogen bonds between ASN318, SER103 and the antagonist still play important role. It is worth noting that the hydrogen bond of ASN318 with the amide oxygen atom (ASN318-NH • • • O=) is

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

13

Fig. 13. Docking results and 2D non-bonded interaction maps of designed compounds DS01 in the binding site of OX1 R pocket (A, B) and OX2 R pocket (C, D).

Fig. 14. Root-mean-square deviation (RMSD) of the 4ZJ8-33 (marked in blue) and 4ZJ8-DS01 (marked in red) complexes versus the dynamics simulation time. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

˚ which is more stable than the 2.018 A˚ before the dynamic 1.838 A, simulation. The hydrogen bond formed by the nitrogen atom at the 2-position of the pyrimidine ring and residue ASN318 disappeared. The nitrogen atom at 6-position of pyrimidine ring formed a new 2.194 A˚ hydrogen bond with residue GLN126 (GLN126-NH • • • N). Therefore, the dynamics simulation results demonstrate the importance of hydrogen bonds for the antagonists.

3.6. ADME analysis ADME prediction plays a significant role in the process for new drugs since the physicochemical properties of the drug can directly affect the pharmacokinetic results. Prediction statistical parameters results of the new design compounds DS01 and DS02 via the online network tool SwissADME were illustrated in Table 5. The

14

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

Fig. 15. Computational modeling of 4ZJ8-33 (A) and 4ZJ8-DS01 (B) binding. (A) Overlay of initial structure (green) and final structure (orange) from 10-ns simulation in OX1 R pocket. (B) Overlay of initial structure (mangeta) and final structure (orange) from 10-ns simulation in OX1 R pocket. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 5 Structures of designed molecules and their predicted pIC50 based on the OX1 R and OX2 R models. NO.

MW (g/mol)

LogP

TPSA (A˚ 2 )

Log S

Num. rotatable bonds

HIA

BBB

CYP2C9 inhibitor

Synthetic accessibility

Optimal scope 33 DS01 DS02

150–500 449.46 467.45 492.48

−0.7–5.0 3.60 3.42 2.93

20–130 89.69 89.69 114.10

≤6 −5.01 −5.17 −5.29

≤9 8 8 9

– High High High

– No No No

– Yes Yes Yes

1–10 (easy-difficult) 3.81 3.90 4.19

logP values of the designed representative compounds DS01 and DS02 were 3.42 and 2.93, respectively. The logS values were −5.17 and −5.29, respectively. It indicated that the new selective OX1 R antagonists not only has good absorption capacity, but also possessed reasonable solubility properties in vivo. The TPSA (Topological polar surface area) values were 89.69 A˚ 2 and 114.10 A˚ 2 , respectively, indicating that DS01 and DS02 had excellent transport properties in vivo. Two compounds all exhibit high HIA (Human gastrointestinal absorption) and the yolk (yellow region) represents the high probability of brain penetration, which was the same as compound33. The “Yes” in Table 5 indicated that the compound has a greater probability of being the inhibitor of CYP2C9. Therefore, the newly designed compounds can be excreted via metabolic biotransformation by the inhibition of the cytochrome CYP2C9 enzyme. In addition, it is noteworthy that the synthetic accessibility of the designed compounds DS01 and DS02 were 3.90 and 4.19, respectively, indicating that they were not difficult to synthesize. In short, pharmacokinetics predictions could improve the success rate of newly designed antagonists.

Table 6 QFIT values and molecular docking scores of Suvorexant, molecular 33, and screened compounds. QFIT

Total scorea (4zj8)

Suvorexant



5.5668

33



7.6619

ZINC77124457

42.670

9.0412

ZINC01751692

46.290

8.9842

ZINC01364054

42.910

8.4011

ZINC09711923

45.080

7.8287

NO.

Structure

3.7. Discovery of new OX1 R antagonists by virtual screening A series of pharmacodynamic elements of compounds are called a pharmacophore. The pharmacophore information cannot only provide structural features of the ligand, but also can be used to screen new skeleton drug molecules. The optimal pharmacophore model (Model-001, Specificity = 5.620) was presented in Fig. 16. There were eight pharmacophores features, including four green pharmacophore features, which indicate the role of hydrogen bond receptors. Two of them were located in the R1 region, the third was located in amide oxygen atom at linkage moiety and the last one was in the R4 region. The remaining four blue–green were hydrophobic interaction pharmacophore features which located in the four rings that form the molecular skeleton of the pharmacophore. The four compounds in Table 6 (ZINC77124457, ZINC01751692, ZINC01364054, and ZINC09711923) were screened from the ZINC12 database of 1020,679 compounds, which the QFIT values were all over 40. The Suvorexant scored 5.5668 and 6.3659 for the OX1 R and OX2 R, respectively, and the compound 33 scores were 7.6619

a

Total score: binding affinities of protein-ligand complexes.

and 7.1655, respectively. The four compounds we screened scored 9.0412, 8.9842, 8.40111 and 7.8287, respectively, which were higher than the original ligands (Suvorexant) and compound 33. And they scored 5.8777, 7.7287, 5.0842 and 7.0058 in OX2 R, which were lower than the score in OX1 R. Docking results hydrogen bond pockets of the four new skeleton compounds by virtual screen were

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

15

Fig. 16. Pharmacophore model features for virtual screening. Three were four hydrogen bond acceptors sites (AA-1, AA-2, AA-3 and AA-4, colored in blue), four hydrophobic interaction sites (HY-5, HY-6, HY-7 and HY-8, colored blue–green). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 17. Four new skeleton compounds based on virtual screening (ZINC77124457, ZINC01751692, ZINC01364054 and ZINC09711923). A, B, C and D show docking results of four molecular into the hydrogen bond active site of OX1 R, respectively. E, F, G and H indicate docking results of four molecular into the hydrogen bond active site of OX2 R, respectively.

shown in Fig. 17. It was found that the molecules (ZINC77124457, ZINC01751692, ZINC01364054, and ZINC09711923) not only form a significant hydrogen bond with residue ASN318, but also form additional hydrogen bonds, such as SER103, GLN179, ARG322, CYS202, and HIS216. In OX2 R, hydrogen bonds were formed with residue ASN324, ARG328, MET191, and LYS327. From Jie Yin’s com-

putational analysis, the interaction of the SER103 with the antagonist in the OX1 R may play an important role for the selectivity for OX1 R over OX2 R. It is worth noting that the screened compounds (ZINC77124457 and ZINC09711923) formed hydrogen bonds with SER103 in the OX1 R. Hence, the four screened compounds may be potential selective OX1 R antagonists.

16

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17

4. Conclusion In this study, a systematic theoretical structure-activity relationship of 36 pyrazoylethylbenzamide derivatives was discussed by computer-assisted drug method, including QSAR models, molecular docking, molecular dynamics, ADME prediction, pharmacophore model, and virtual screening. OX1 R QSAR models were established (CoMFA: q2 = 0.788, r2 = 0.988, r2 pred = 0.787; CoMSIA: q2 = 0.783, r2 = 0.953, r2 pred = 0.817; HQSAR: q2 = 0.512, r2 = 0.927, r2 pred = 0.917), and the results indicated that three constructed models were reliable and can be used to predict the biological activity of newly designed compounds. Contour maps revealed that hydrogen bond receptor was one of the main factors affecting the biological activity of antagonists in the OX1 R models, especially the nitrogen atom on pyrimidine ring 2 position, the amide oxygen atom, and the benzene ring-linked fluorine atom. Furthermore, the stable binding of the antagonist to the pocket by molecular docking can be attributed to hydrogen bonds between the antagonists and residue ASN318 and SER103. To design more selective OX1 R antagonists, we introduced a fluorine atom or a carboxyl group at the 6 position of the fluorobenzene ring in the R1 region and a hydroxyl group at the 2 position of benzene ring in the R4 region of the template molecule 33. Finally, seven new 1-SORAs (DS01–DS07) were designed. The results of ADME prediction showed that the pharmacokinetic parameters (LogP, Log S, HIA, CYP2C9 inhibitor) of the designed compounds not only meet the required pharmacokinetic parameters but also display potential druggability properties and higher selectivity index than the template compound 33 (SI = 265). Hydrogen bonds between the nitrogen atom in the R1 region, the fluorine atom in the R4 region and the OX1 R pocket were confirmed by molecular dynamic simulations. Furthermore, OX2 R QSAR models (CoMFA: q2 = 0.852, r2 = 0.979, r2 pred = 0.846; CoMSIA: q2 = 0.823, r2 = 0.965, r2 pred = 0.797; HQSAR: q2 = 0.654, r2 = 0.951, r2 pred =0.894) were constructed and compared with the OX1 R QSAR models. In addition, to discover more novel OX1 R antagonist, four potential new skeletons antagonists (Z1–Z4) were screened from the ZINC12 database containing 1020,679 compounds based on the best pharmacophore model (Model-001, Specificity = 5.620). The screened new backbone compounds (1,2,4-oxadiazole and quinoxaline skeleton) formed important hydrogen bond with OX1 R residues (ASN318 and SER103) and displayed additional hydrogen bond (GLN176 and ARG322) to stabilize the compounds. These theoretical results may provide guidance for the design, optimization, and synthesis of more selective OX1 R antagonists. New designed compounds (DS01–DS07) and screened molecules (Z1–Z4) may be candidate compounds for the relevant research groups. The results might also provide an important strategy for developing insomnia compounds. Acknowledgments We acknowledge the Shanghai Science and Technology Commission for support of this work. Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.jtice.2019.03.018. References [1] Buysse DJ. Insomnia. JAMA J Am Med Assoc 2013;309:706–16. [2] Ohayon MM. Epidemiology of insomnia: what we know and what we still need to learn. Sleep Med Rev 2002;6:97–111.

[3] Taylor DJ, Mallory LJ, Lichstein KL, H Heith D, Riedel BW, Bush AJ. Comorbidity of chronic insomnia with medical problems. Sleep 2007;30:213–18. [4] Sibyl A, Inge P, Hilde H, Pascale S, Tom D, Thierry C. Barriers to nonpharmacologic treatments for stress, anxiety, and insomnia: family physicians’ attitudes toward benzodiazepine prescribing. Can Fam Physician 2010;56:398–406. [5] Holbrook AM, Crowther R, Lotter A, Cheng C, King D. Meta-analysis of benzodiazepine use in the treatment of insomnia. J Am Geriatr Soc 2010;49:824–6. [6] Wagner J, Wagner ML. Non-benzodiazepines for the treatment of insomnia. Sleep Med Rev 20 0 0;4:551–81. [7] Huedo-Medina TB, Kirsch I, Middlemass J, Klonizakis M, Siriwardena AN. Effectiveness of non-benzodiazepine hypnotics in treatment of adult insomnia: meta-analysis of data submitted to the Food and Drug Administration. BMJ 2013;345:8343–55. [8] Saper CB, Scammell TE, Jun L. Hypothalamic regulation of sleep and circadian rhythms. Nature 2005;437:1257–63. [9] Morin CM, Drake CL, Harvey AG, Krystal AD, Manber R, Riemann D, et al. Insomnia disorder. Nat Rev Dis Prim 2015;1:15026–43. [10] de Lecea L, Kilduff T, Peyron C, Gao X-B, Foye PE, Danielson PE, et al. The hypocretins: hypothalamus-specific peptides with neuroexcitatory activity. P Natl Acad Sci USA 1998;95:322–7. [11] Sakurai T, Amemiya A, Ishii M, Matsuzaki I, Chemelli R, Tanaka H, et al. Orexins and orexin receptors: a family of hypothalamic neuropeptides and g protein– coupled receptors that regulate feeding behavior. Cell 1998;92:573–85. [12] Smart D, Brough SJ, Jewitt F, Johns A, Porter RA, Jerman JC. SB-334867-A: the first selective orexin-1 receptor antagonist. Brit J Pharmacol 2001; 132:1179–82. [13] Roecker AJ, Cox CD, Coleman PJ. Orexin receptor antagonists: new therapeutic agents for the treatment of insomnia. J Med Chem 2016;59:504–30. [14] Roecker A, Coleman P, Roecker AJ, Coleman PJ. Orexin receptor antagonists: medicinal chemistry and therapeutic potential. Curr Top Med Chem 2008;8:977–87. [15] Kuriyama A, Tabata H. Suvorexant for the treatment of primary insomnia: a systematic review and meta-analysis. Sleep Med Rev 2017;35:1–7. [16] Rhyne DN, Anderson SL. Suvorexant in insomnia: efficacy, safety and place in therapy. Ther Adv Drug Saf 2015;6:189–95. [17] Black S, Morairty S, Fisher S, Chen T-M, R Warrier D, Kilduff T. Almorexant promotes sleep and exacerbates cataplexy in a murine model of narcolepsy. Sleep 2013;36:325–36. [18] Wang C, Moseley C, M Carlin S, Wilson C, Neelamegam R, Hooker J. Radiosynthesis and evaluation of [C]EMPA as a potential PET tracer for orexin 2 receptors. Bioorg Med Chem Lett 2013;23:3389–92. [19] Rodgers RJ, Halford JC, Rl NDS, Al CDS, Piper DC, Arch JR, et al. SB-334867, a selective orexin-1 receptor antagonist, enhances behavioural satiety and blocks the hyperphagic effect of orexin-A in rats. Eur J Neurosci 2010;13:1444–52. [20] Steiner MA, Gatfield J, Brisbare-Roch C, Dietrich H, Treiber A, Jenck F, et al. Discovery and characterization of ACT-335827, an orally available, brain penetrant orexin receptor type 1 selective antagonist. ChemMedChem 2013;8:898–903. [21] Futamura A, Nozawa D, Araki Y, Tamura Y, Tokura S, Kawamoto H, et al. Identification of highly selective and potent orexin receptor 1 antagonists derived from a dual orexin receptor 1/2 antagonist based on the structural framework of pyrazoylethylbenzamide. Bioorganic Med Chem 2017;25:5203–15. [22] Patel S, Patel B, Bhatt H. 3D-QSAR studies on 5-hydroxy-6-oxo-1,6-dihydropyrimidine-4- carboxamide derivatives as HIV-1 integrase inhibitors. J Taiwan Inst Chem Eng 2016;59:61–8. ˙ zy ˙ nska-Granica ´ [23] Zy B, Trzaskowski B, Niewieczerzał S, Filipek S, Zegrocka-Stendel O, Dutkiewicz M, et al. Pharmacophore guided discovery of small-molecule interleukin 15 inhibitors. Eur J Med Chem 2017;136:543–7. [24] Floresta G, Rescifina A, Marrazzo A, Dichiara M, Pistarà V, Pittalà V, et al. Hyphenated 3D-QSAR statistical model-scaffold hopping analysis for the identification of potentially potent and selective sigma-2 receptor ligands. Eur J Med Chem 2017;139:884–91. [25] Jaiteh M, Zeifman A, Saarinen M, Svenningsson P, Brea JM, Loza MI, et al. Docking screens for dual inhibitors of disparate drug targets for Parkinson’s disease. J Med Chem 2018;2017:5269–78. [26] Golbraikh A, Tropsha A. Predictive QSAR modeling based on diversity sampling of experimental datasets for the training and test set selection. J Comput Aid Mol Des 2002;5:357–69. [27] Yin J, Babaoglu K, Brautigam CA, Clark L, Shao Z, Scheuermann TH, et al. Structure and ligand-binding mechanism of the human OX1 and OX2 orexin receptors. Nat Struct Mol Biol 2016;23:293–301. [28] Jie Y, Juan Carlos M, Peter K, Rosenbaum DM. Crystal structure of the human OX2 orexin receptor bound to the insomnia drug suvorexant. Nature 2014;519:247–50. [29] Wu S, Qi W, Su R, Li T, Lu D, He Z. CoMFA and CoMSIA analysis of ACE-inhibitory, antimicrobial and bitter-tasting peptides. Eur J Med Chem 2014;84:100–6. [30] Wang A, Yang Y, Jun Y, Wang B, Lv K, Liu M, et al. Synthesis, evaluation and CoMFA/CoMSIA study of nitrofuranyl methyl N-heterocycles as novel antitubercular agents. Bioorganic Med Chem 2018;26:2073–84. [31] Hirdesh K, Rajendra K, Grewal BK, Sobhia ME. Insights into the structural requirements of PKCβ II inhibitors based on HQSAR and CoMSIA analyses. Chem Biol Drug Des 2011;78:283–8. [32] Seel M, Turner DB, Willett P. Effect of parameter variations on the effectiveness of HQSAR analyses. Qsar Comb Sci 2015;18:245–52.

Q. Li, C. Zhang and Y. Ren / Journal of the Taiwan Institute of Chemical Engineers 100 (2019) 1–17 [33] Tenenhaus M, Vinzi VE. PLS regression, PLS path modeling and generalized Procrustean analysis: a combined approach for multiblock analysis. J Chemom 2010;19:145–53. [34] Alexander DLJ, Tropsha A, Winkler DA. Beware of R2: simple, unambiguous assessment of the prediction accuracy of QSAR and QSPR models. J Chem Inf Model 2015;55:1316–22. [35] Roy P.P., Paul S., Mitra I., Roy K. On Two novel parameters for validation of predictive QSAR models 2009;14:1660–1701. [36] Mitra I, Roy PP, Kar S, Ojha PK, Roy K. On further application of r2m as a metric for validation of QSAR models. J Chemom 2010;24:22–33. [37] Ferreira LG, Santos Dos RN, Glaucius O, Andricopulo AD. Molecular docking and structure-based drug design strategies. Molecules 2015;20:13384–421. [38] Spitzer R, Jain AN. Surflex-dock: docking benchmarks and real-world application. J Comput Aided Mol Des 2012;26:687–99. [39] Rose PW, Prlic´ A, Altunkaya A, Bi C, Bradley AR, Christie CH, et al. The RCSB protein data bank: integrative view of protein, gene and 3D structural information. Nucleic Acids Res 2017;45:271–81. [40] Case DA, Cheatham TE, Tom D, Holger G, Ray L, Merz KM, et al. The Amber biomolecular simulation programs. J Comput Chem 2010;26:1668–88.

17

[41] Ohtsuki S, Uchida Y, Kubo Y, Terasaki T. Quantitative targeted absolute proteomics-based adme research as a new path to drug discovery and development: methodology, advantages, strategy, and prospects. J Pharm Sci 2011;100:3547–59. [42] Daina A, Michielin O, Zoete V. SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci RepUK 2017;7:42717–29. [43] Cheng T, Li Q, Zhou Z, Wang Y, Bryant SH. Structure-based virtual screening for drug discovery: a problem-centric review. AAPS J 2012;14:133–41. [44] Chiang Y-K, Kuo C-C, Wu Y-S, Chen C-T, Coumar MS, Wu J-S, et al. Generation of ligand-based pharmacophore model and virtual screening for identification of novel tubulin inhibitors with potent anticancer activity. J Med Chem 2009;52:4221–33. [45] Irwin JJ, Teague S, Mysinger MM, Bolstad ES, Coleman RG. ZINC: a free tool to discover chemistry for biology. J Chem Inf Model 2012;52:1757–68. [46] Cox CD, Mcgaughey GB, Bogusky MJ, Whitman DB, Ball RG, Winrow CJ, et al. Conformational analysis of N,N-disubstituted-1,4-diazepane orexin receptor antagonists and implications for receptor binding. Bioorg Med Chem Lett 2009;19:2997–3001.