Discovery of novel insomnia leads from screening traditional Chinese medicine database

Insomnia is a prominent modern disease that affects an increasing population. Undesirable side effects of commercial drugs highlight the need to develop novel insomnia drugs. Virtual screening of traditional chinese medicine Database@Taiwan (TCM Database@Taiwan) identified 2-O-Caffeoyl tartaric acid (1), 2-O-Feruloyl tartaric acid (2), and Mumefural (3) as potential agonists for both gamma-amino butyric acid (GABA) or benzodiazepine (BZ) binding sites. The TCM candidates exhibited higher affinity than GABA and Zolpidem, a phenomenon that could be attributed to higher quantity of stabilizing H-bonds. Efficacy profiles using support vector machines and pharmacophore contour also suggest drug potential of the TCM candidates. Fragments added to the de novo derivatives 3a, 3b, 3c for GABA binding site, and 1a, 2a, and 3d for BZ binding site contributed to new binding sites and structural stability, further optimizing binding to GABA or BZ binding sites. Increased opening of the ion channel by candidate ligands provide strong support for their potential biological functions. The dual binding properties of the TCM candidates present a unique opportunity to develop twin-targeting drugs with less side effects. Derivative structures can be used as starting points for developing high affinity GABAA receptor agonists with specificity towards GABA binding site and BZ binding site.


Introduction
Insomnia is a serious modern illness affecting approximately 10-38% of the population (Dundar et al., 2004;Whiting, 2006). Hyperarousal of the central nervous system (CNS) is associated with insomnia. Gamma-amino butyric acid (GABA) is an important substance in the cerebral cortex for countering neural overexcitation and inhibiting neural signal transmission (Alvarez Dolado & Broccoli, 2011). The restorative effect of peaceful sleep is linked to GABA production.
GABA generally regulates neural excitation by binding to either GABA A receptor (GABA A R) or GABA B receptor (GABA B R), and eliciting changes in ion concentrations. GABA A R is associated with the increase of Cl À ions (Olsen & Sieghart, 2009), and GABA B R activates K + channels and inhibits Ca 2+ channels (Ulrich & Bettler, 2007). With regard to sleep and insomnia treatment, activation of GABA A R by GABA will open the Cl À ion channel, leading to a membrane potential change by the increased Cl À concentration and inhibit neural hyper-polarization (Rudolph & Knoflach, 2011). GABA A R are penta-hetero membrane proteins formed by 19 identified subunits (α1-α6, β1-β3, γ1-γ3, δ, ɛ, π, h, and ρ1-ρ3) (Kuriyama, Hirouchi, & Nakayasu, 1993). The most common conformation found in CNS is α1-β2-α1-γ2-β2. Two key binding sites, GABA binding site and benzodiazepine (BZ) binding site, for regulating neural hyperarousal have been identified. GABA binding site is located at the interface of α1/β2 on GABA A R, and BZ binding site is located at the interface between α1/γ2 on GABA A R (Sarto-Jackson & Sieghart, 2008;. Many commercial sedatives or sleep medications contain GABA-like active ingredients to adjust for insomnia which occurs under low GABA concentrations. These hypnotics can be classified into benzodiazepines (BZD) and non-benzodiazepines (NBZ) according to structural characteristics (Nutt & Stahl, 2010;Rosenberg, 2006). Commonly used treatments for insomnia are listed in Table 1. Though successful to varying degrees for relieving insomnia, side effects are commonly reported and can be severe (Bocca & Denise, 2000;Kales, Scharf, Kales, & Soldatos, 1979;Kirkwood, 1999;Kovacic & Somanathan, 2009;Lader, 2001;Monchesky, Billings, & Phillips, 1986;Patat, Paty, & Hindmarch, 2001;Pierce, Shu, & Groves, 1990), resulting in the demand for safer insomnia relief. Recently, many researches attempt to detect some agents acted on GABA A R for relieving insomnia without side effects Moree et al., 2010;Winrow et al., 2012).
Previous work in our lab shows the potential of traditional Chinese medicine (TCM) as an important source for new treatment against insomnia (Chen, 2007(Chen, , 2008(Chen, , 2009Chen, Chen, Wu, & Tsai, 2008). From the experimental scheme outlined in Figure 1, we aimed to investigate ligands with dual targeting properties by utilizing a double fold virtual screening process for GABA binding site and BZ binding site. De novo evolution was employed to derive chemical structures for optimum binding affinities.

Homology modeling of GABA A receptor
Homology modeling was performed to construct a model structure due to the lack of GABA A R crystal structures in Protein Data Bank. Subunit sequences downloaded from Swiss-Prot were P14867 for α1, P47870 for β2, and P18507 for γ2. Sequences were arranged to form the major combination type α1β2α1γ2β2 (Rudolph & Knoflach, 2011). Amino acid residues extracted from each sequence for modeling were 28-251 in α1, 25-244 in β2, and 40-273 in γ2 (Chou, Watenpaugh, & Heinrikson, 1999). The pentamer acetylcholine binding protein (AChBP) from Lymnaea stagnalis (PDB: 1I9B) (Brejc et al., 2001) was used as the reference structure. GABA A R structure was modeled using the Build Homology Models protocol in Discovery Studio 2.5 (DS 2.5) BZD Benzodiazepine withdrawl syndrome could result from discontinuation. Symptoms include seizures, psychosis, and severe forms of anxiety and insomnia (Kales et al., 1979). Zaleplon (Sonata ® and Starnoc ® ) NBZ Dizziness and headaches (Kirkwood, 1999;Lader, 2001;Patat et al., 2001).

Virtual screening and de novo evolution
The modeled GABA A R structure was used to virtually screen for ligands from TCM Database@Taiwan (http:// tcm.cmu.edu.tw/) (Chen, 2011), that could bind to GABA and BZ binding sites on GABA A R. All ligands used for screening passed Lipinski's Rule of Five. GABA and BZ binding sites were defined as the space between subunits α1 and β2 (GABA binding site), and α1 and γ2 (BZ binding site), respectively Klausberger, Fuchs, Mayer, Ehya, & Sieghart, 2000;. Docking was performed with LigandFit under the forcefield of Chemistry at HARvard Molecular Mechanics (CHARMm) (Brooks et al., 2009). GABA and Zolpidem were adopted as the respective controls for GABA binding site and BZ binding site. Twenty TCM ligands with the highest Dock Score for each binding site were selected for de novo evolution. Novel derivatives were generated by docking the TCM ligands into their respective receptor site and adding new fragments from the fragment library in DS 2.5 based on the chemical properties of the binding site. Full evolution was selected with a binding site sphere of 7 . Derivatives were ranked according to Dock Scores. The top three TCM ligands and derivatives for each binding site were selected as candidates for further analysis.

Bioactivity prediction using support vector machines
Due to the lack of available structure-bioactivity correlation research for GABA binding site, quantitative structure activity relationship (QSAR) prediction model was only constructed for BZ binding site. Structures and biological affinities of 26 BZ binding site agonists (Grunwald et al., 2006) were subjected to Genetic Function Approximation (GFA) (Yang, Chang, Chen, & Chen, 2011) in DS 2.5 to determine representative molecular descriptors associated with bioactivity. Representative descriptors were used to construct the support vector machines (SVM) model through LibSVM program (Yang et al., 2010). Models were subjected to fivefold cross validation and strength of the generated prediction models were evaluated by square correlation coefficient (R 2 ). The model with the highest cross validation and R 2 values was used to predict biological affinities of the selected candidates for BZ binding site.

Molecular dynamics simulations
To ascertain stability of candidates in GABA binding site and BZ binding site, molecular dynamics (MD) of each protein-ligand complex was conducted. The MD simulations protocol was conducted using Standard Dynamics Cascade module and Dynamics (Production) module in the Simulation package of DS 2.5. The protocol of the minimization steps were conducted by both Steepest Descent (Fletcher, 1969) and Conjugate Gradient (Fletcher & Reeves, 1964) with max number of 500 cycles, respectively. Then, the heating step gradually heated the whole system from initial temperature of 50 K to target temperature of 310 K over 20 ps. After a following 100 ps equilibration step, a 20 ns MD was conducted with fixed H-atoms under NVT canonical ensemble using CHARMm forcefield (Brooks et al., 2009) and time steps of 1 fs. For the Berendsen thermal coupling method, the temperature coupling decay time was set to 5 ps. A total of 2000 snapshots were taken at 10 ps intervals of each protein-ligand complex. Snapshots were used to analyze root mean square deviations (RMSD), total energy, H-bond distance/quantity, and torsion fluctuations. Ligand pathways were calculated using the LigandPath program developed by Dr Tu-Liang Lin (Lin & Song, 2011). The free space calculated for conformations from 0 to 5 ns were defined as "entry pathways", and that for conformations from 16 to 20 ns were defined as "exit pathways". Minimum clearance and surface probe radius were set at 3 and 5 Å, respectively. Channel opening was used to visualize the functional effect of candidate binding to GABA binding site and/or BZ binding site.

Pharmacophore features
To further provide drug development insight, hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), and hydrophobic (HY) features were calculated for GABA binding site and BZ binding site. The site sphere radius was set at 10 Å for GABA binding site and 12 Å for BZ binding site. The selected candidates were docked into the respective binding sites to observe ligand contour with the pharmacophore features of the binding site.

Results and discussion
Homology modeling of GABA A R Sequence alignment between AChBP and GABA A R show a sequence identity of 19.8% and a sequence similarity of 45.1% (Supplementary Figure S1). A comparison between reported key residues in the literature and the corresponding residues in this study are given in Supplementary Table S1. Spatial location of the GABA binding site and BZ binding site are illustrated in Fig Figure 3. The top three compounds from TCM ligands and derivatives per binding site selected for further analysis in addition to the control compounds are summarized in Figure 3. For clarification purposes, ligands screened directly from the TCM Database@Taiwan are termed "TCM candidates", and derived compounds are referred to as "derivatives". Derivatives are named according to the number designated for their parent TCM compound followed by a lower case alphabet. TCM candidates and derivatives are collectively referred to as candidate ligands. Interestingly, 2-O-Caffeoyl tartaric acid (1), 2-O-Feruloyl tartaric acid (2), and Mumefural (3) exhibited the highest Dock Scores in both GABA and BZ binding site. This suggests a unique possibility for designing single drugs with dual targets. As expected, derivatives exhibited higher docking score properties (Dock Score, -PLP1, -PLP2, -PMF) than their parent TCM candidate, implying increased receptor affinity through the added structural fragments. The higher affinity of candidate ligands compared to GABA and Zolpidem may be due to electrostatic and van der Waals (vdW) interactions.
Receptor-ligand interactions during docking are summarized in Table 2. GABA was stabilized by H-bonds with Glu365 and Lys408, vDW interactions with Leu311, Val410, Phe412, and Val613, and HY interaction with Phe412. By comparison, candidate ligands formed more interactions with the GABA binding site. Tyr309, Arg419, Arg499, Arg517, Arg552, and Arg605 formed H-bonds with all candidate ligands, suggesting importance for ligand binding. In addition, Tyr369, Tyr417, Leu560, and Thr562 are also key residues which form HY contact with at least five of the candidate ligands. De novo evolution of 3 enabled additional interactions with Phe412, Phe478, Phe497, Thr562, and Val613, all of which might have contributed to the increased Dock Score observed for 3a, 3b, and 3c in Table 3. In BZ binding site, Zolpidem was stabilized by a single H-bond with Arg852 while candidate ligands were anchored by multiple H-bonds (Table 2). H-bond interactions with Arg799 and Arg852 were observed in all candidate ligands. HY interactions between Zolpidem and the candidate ligands were not significantly different; nonetheless, increased HY interactions were still observed when comparing 1a to 1 and 3d to 3. His534, Val635, and Tyr642 appear to be important HY residues in BZ binding site. Arginine was the primary type of amino acid candidate ligands interacted with. The carbonyl moieties located on candidate ligands could have higher affinity to arginine, leading to the formation of Hbonds and electrostatic interactions. Effect of the carbonyl groups is also reflected directly on DockScore. The higher DockScores of candidate ligands indicate better electrostatic interactions and vdW with the receptor, which could be attributed to the increased binding ability caused by the carbonyl groups. Docking results indicate that ligand candidates exhibit desirable affinity qualities based on their ability to form more stabilizing interactions than either GABA or Zolpidem in their respective binding sites.

Bioactivity prediction using SVM
Representative molecular descriptors determined by GA for BZ binding site agonists include Num_Chain Assemblies, Molecular_PolarSurfaceArea, IAC_Mean, Kappa_2_AM, Dipole_Z, Jurs_WPSA_3, and Shadow_-XYfrac. Definitions of each descriptor indicate that bioactivity was determined to be linked primarily to the size polarity, and shape of the ligand (Supplementary Table S6). The correlation between observed and predicted pK i of R 2 = 0.7313 indicates that the constructed SVM model is of acceptable accuracy (Supplementary Figure S2). Predicted biological affinities (pK i ; given in parenthesis) for BZ binding site candidate ligands, from high to low are: 1 (7.21) > 1a (7.20) > 3 (7.20) > 3d (7.19) > 2a (7.18) > Zolpidem (6.27) > 2 (4.53). Candidate ligands for BZ binding site were predicted to have higher biological affinities with the exception of 2-O-Feruloyl tartaric acid (2).

MD simulations
Trajectories of complex RMSD and total energy indicate that equilibrium was achieved for all complexes by the end of MD (Figure 4). Total energies of equilibrated receptor-ligand complexes in GABA binding site ranged between À58,000 and À59,000 kcal/mol. Though no clear pattern could be established among TCM candidates and derivatives with regard to total energy or complex RMSD, candidate ligands with higher ligand and complex RMSDs stabilized at lower total energies. Higher RMSD values could reflect more significant changes which lead to lower energy conformations. In BZ binding site, complexes formed by ligand candidates had lower total energies compared to the Zolpidem complex. No clear relationship between RMSD and total energies could be formed, but all derivatives complexes had lower total energies than the TCM complexes.
Sustained presence of an H-bond during MD can be viewed as an indicator of H-bond stability. Frequencies of H-bonds formed by candidate compounds with GABA binding site and BZ binding site are shown in Tables 4 and 5, respectively. To reduce under estimation of occupancies due to the occupancy cutoff of 2.5 Å, low occupancy H-bonds were further verified with H-bond trajectories.
In GABA binding site, 1 formed high occupancy Hbonds with Tyr417, Arg499, Arg517, and Arg552 (Table 4). A longer distance (hence lower occupancy) H-bond was maintained with Arg605 as well ( Figure 5). The H-bond with Arg419 was lost due to increased distance from 3.0 to 4.5 Å. 2 was stabilized by H-bonds with Tyr309, Arg499, Arg517, and Arg552 during MD. The major difference from docking (Table 2) was that Hbonds with Tyr417, Arg419 or Arg605 were not observed. 3 was anchored in GABA binding site by all H-bonds observed during docking except for Tyr309. Hbonds with Arg419 and Arg605 were weaker due to the longer H-bond distance ( Figure 5). Compared to docking, an additional H-bond with Tyr417 was observed. MD reveals that residues Arg499 and Arg552 are critical residues that could form stabilizing H-bonds with TCM candidates in GABA binding site. Two important conclusions can be made from derivative results. Since all derivatives formed stable interactions with Arg499 and Arg552, importance of these two residues for ligand binding with GABA binding site is further highlighted. Secondly, added fragments on 3a, 3b, and 3c increased the amount of stabilizing interactions formed to anchor the derivatives. In addition to H-bonds observed with the native Mumefural structure (Tyr417, Arg419, Arg499, Arg552, Arg605), additional stable H-bonds were formed with Tyr309, Thr561, and Thr562 in 3a, Tyr309, Arg517, and Thr562 in 3b, and Tyr309 in 3c ( Figure 5).
Likewise, key binding residues for BZ binding site can be investigated. 2 formed stable interaction with His534, Gly536, Arg799, and Arg852 during MD (Table 5). H-bond with Lys588 was not stable (Figure 6). 1 maintained H-bonds with Arg799 and Arg852, and formed an additional H-bond with His534. Interaction with Tyr592 or Ser637 was not detected during MD. 3    Figure 7 also provide important information regarding H-bond stability. Among the torsions between 1 and GABA-binding site that are related to H-bond formation, torsion (2) is the most unstable, and this unstable torsion is probably the underlying reason for the inability of 1 to form stable H-bond with Arg419 (Table 4). No significant changes were observed for (5), therefore H-bonds with Arg222 and Arg517 are very stable. In 2, large fluctuations at (8) may have caused low occupancy with Tyr396, whereas the relatively stable (10) enables high occupancy H-bonds with  Tyr417, Arg499, and Arg517. Torsion (12) in 3 affects stable interaction with Arg419. Stabilizing of (16), after 5 ns, enabled a longer distance H-bond with Arg605 ( Figure 5). Locations on 3a associated with H-bonds are (17) and (21). Torsion fluctuations at these locations are small, forming stable H-bonds with Tyr309, Tyr417, Arg419, Arg517, and Arg522. The added fragment of 3a (at torsion (17)) forms a high occupancy (99.95%) H-bond with Arg419. leading to higher stability of the derivative compared to 3. Torsion trajectories at (23) in 3b exhibits similar trends to H-bond distance fluctuations at Tyr309. The 67% H-bond occupancy could be  attributed to the large fluctuations of (23) during the first 10 ns, which hinder the successful formation of H-bonds. Small variations at (24) and (25) enable stable interactions to be formed with Tyr417, Arg517, Arg552, and Thr562 during MD. Torsion changes in 3b were generally smaller than those for 3. The higher stability could be related to additional H-bonds formed with Tyr309 and Thr569. H-bond related torsion locations measured for 3c are identical to those for Mumefural. The stable, near identical fluctuation trajectories between 3c and 3 could provide some explanation for the similar complex total energy observed in Figure 4. In BZ binding site, instability of (34) affected Hbond formation of 2 with Lys588. (36) and (37) showed little variation, therefore the H-bonds were stable. Similarly, (41) in 1 was stable, forming a high occupancy Hbond with Arg799; (42) fluctuated, leading to large distance changes to Arg852 and a lower H-bond occupancy. Despite fluctuations at (47) in 3, H-bonds with Arg799 and Arg852 were unaffected. In 1a, torsions related to  H-bonds ((48), (51), (52), (53)) were stable, implying that 1a was well anchored within the complex. The lower occupancy (14.05%) ( Table 5) and larger H-bond variation observed between 2b and Lys588 ( Figure 6) might be related to variations at (54). Torsions (63) and (64) in 3d are associated with H-bonds at Gly536, Ser637, Arg799, and Arg852. Since these two torsions were very stable, H-bonds with the aforementioned residues were also stable throughout MD. Limited fluctuations at (67) and (68) were observed for Zolpidem, suggesting a stable H-bond with Gly536 throughout MD.
Torsions associated with added fragments can also provide critical information on the stabilizing effect of said fragments (Figure 7). Within GABA binding site, the HY aromatic ring of 3 is located at (13) with fluctuations of 70°. Corresponding angles at 3a, 3b, and 3c are located at (18), (24), and (29), respectively. The torsion changes at these locations are smaller than that at (13), suggesting that the new moieties, added by de novo evolution, add stability to the derivative-GABA binding site complex. For BZ binding site, the HY aromatic rings for 2, 1, and 3 are located at (35), (40), and (44), with corresponding fluctuations of 80°, 85°, and 105°. 1a, 2a, and 3d at (51), (56), and (61) show significantly smaller fluctuations than (35), (40), and (44). Results indicate that HY moieties in the de novo derivatives are more stable than those in their native TCM candidates. In both GABA binding site and BZ binding site, derivatives have more stable binding poses than the native TCM candidates.
As further assessment of drug potential, accessibility of candidate ligands to respective binding sites were evaluated. Derivatives generated more pathways than TCM candidates in both GABA ( Figure 8) and BZ binding sites (Figure 9), with the exception of 2a in BZ binding site (Figure 9(E)), suggesting higher accessibility and potential for forming protein complexes. Zolpidem was found to have three entry and three exit pathways to BZ binding site (Figure 10). TCM candidates had slightly lower accessibility to BZ binding site than Zolpidem (Figure 9(A-C)). However, de novo evolution of 1 and 3  to 1a (Figure 9(D)) and 3d (Figure 9(F)) greatly increased candidate accessibility to BZ binding site.

Effect of structural optimization via de novo evolution on protein-ligand stability
Insights to important binding residues and candidate structures can be extrapolated by cross comparing candidate ligand results in each binding site. Arg517 is a key binding residue for GABA binding site (Figure 11(A)). Derivatives 3a, 3b, and 3c not only maintain original bonds found in 3, but form additional bonds with Tyr309, Arg499, Tyr562, and Arg605. The additional bonds increase stability of derivatives within GABA binding site.
Arg799 and Arg852 are important binding residues for BZ binding site. TCM candidates primarily bound to BZ binding site through the carbonyl groups. As shown in Figure 11(B), due to the increased ring structure in 1a. additional bonds with Arg787 were formed. The added ring structure also increases stability of the native structure, allowing the O atom on the backbone structure to form stable bonds with Arg852. Carbonyl and hydroxyl groups are important binding moieties for 2 and 2a, enabling H-bond formation with Lys588, Arg799, and Lys852 (Figure 11(C)). The added fragment of 3d did not directly provide additional binding sites, but contributed to steric bulk. Nonetheless, H-bonds with Gly536 and Ser637 were observed for 3d in addition to those formed by 3 (Figure 11(D)).

Pharmacophore feature analysis
Pharmacophore features for GABA binding site are shown in Figure 12(A). H-bond related characteristics appear to be primary features. The HBA and HBD pharmacophores detected near Arg499, Arg517, and Arg552 show that this region favors H-bond formation. Our MD results indicate that stable H-bonds were indeed formed between candidate ligands and Arg499, Arg517, and Figure 10. Pathways of Zolpidem to and from BZ binding site during 20 ns MD. Entry (blue) and exit (yellow) pathways are defined by the calculated available space from 0 to 5 and 16-20 ns, respectively. Figure 11. Key ligand candidate binding residues in GABA binding site and BZ binding site. (A) Important binding residues for 3, 3a, 3b, and 3c in GABA binding site. White, blue, orange, and pink represent amino acids that bind to 3, 3a, 3b, and 3c respectively. Binding residues common for all four candidate ligands are shown in green. (B) Important binding residues for 1 and 1a in BZ binding site. Key binding residues for 1 are shown in white, and those for 1a are shown in blue. (C) Important binding residues for 2 and 2a in BZ binding site. Key binding residues for 2 are shown in white, and those for 2a are shown in blue. Common binding residues are shown in green. (D) Important binding residues for 3 and 3d in BZ binding site. Key binding residues for 3 and 3d are shown in white and blue, respectively. Common binding residues are shown in green. White scaffolds and orange fragments represent parent TCM candidate structures and the added fragments of derivative, respectively. Common binding sites for the TCM candidate and its derivative are highlighted with a yellow background.
Arg552. HY features were favored near the Arg419 and Arg499 region. HY moieties of candidate ligands contoured to the HY feature location. In BZ binding site, H-bond related features were detected near His534, Arg799, and Arg852. The imidazo [1,2-a] pyridine of Zolpidem was located near the HBD feature near His534, and the carbonyl group was located near the HBA feature near Arg852. Good contour of Zolpidem with BZ binding site pharmacophores likely contributed to its good observed bioactivity. Superimposition of candidate ligand backbones show that the two structural carbonyls were located near the HBA and HBD pharmacophores. MD results confirm the formation of stable H-bonds with Arg799 and Arg852 by candidate ligands. HY features were located in regions near Leu795. HY moieties of the candidate ligands were found to be adjacent to the HY pharmacophores near Leu795. In summary, contour of candidate ligands to the respective pharamacophores of GABA binding site and BZ binding site imply that the candidate ligands have structural features that are likely to elicit drug-like effects in GABA binding site and BZ binding site.

Effect of candidate binding on ion channel opening
The most important function of GABA A R activation is the opening of ion channels which allows influx of Cl À and neutralization of over-polarization. Though direct experimental results are not attainable in a strictly dry lab environment, channel opening changes observed during MD simulation provides indirect information on biological efficacies of the candidate ligands. Enlargement of the channel opening by TCM candidate binding to GABA binding site is illustrated in Figure 13. No correlation between candidate ranking and intensity of channel opening change was observed. While 1 and 3 generally caused receding of the protein structure along the original opening, a prominent cavity was elicited by  2 (highlighted by yellow square). Derivatives 3a (Figure 13(D)) and 3b (Figure 13(E)) showed additional openings highlighted in the yellow boxes compared to 3. No significant differences were detected between 3 and 3c. When bound to BZ binding site, ion channel changes elicited by the candidate ligands ( Figure 14) were very similar to those by Zolpidem (Figure 15). This implies that the ligand candidates selected for BZ binding site may possess similar functionality to Zolpidem. Comparing between TCM candidates and derivatives, we also speculate that 1a and 2a may have higher efficacy in increasing ion channel opening than 1 and 2, respectively.

Conclusion
Herein we identify TCM compounds, 2-O-Caffeoyl tartaric acid, 2-O-Feruloyl tartaric acid, and Mumefural, with dual binding properties for GABA and BZ binding sites. The TCM candidates exhibited higher affinity towards the two binding sites than GABA and Zolpidem, presenting a unique opportunity to develop dual acting drugs for insomnia. Without the need for multiple drugs that specifically target GABA or BZ binding site, undesirable side effects associated with traditional insomnia drugs could be reduced. We also propose optimized structures 3a, 3b, 3c for GABA binding site, and  structures 1a, 2a, and 3d for BZ binding site. Fragments added to the parent TCM backbone were found to contribute to new binding sites or add to structural stability, leading to better binding characteristics than their parent TCM structure. Efficacy properties including bioactivity prediction by SVM and pharmacophore analysis also suggest drug potential of the candidate ligands. Similar, if not greater, opening of the ion channel by candidate ligands compared to Zolpidem strongly supports the proposed candidates as starting points for the development of new insomnia drugs.