In Silico Screening of Novel α1-GABAA Receptor PAMs towards Schizophrenia Based on Combined Modeling Studies of Imidazo [1,2-a]-Pyridines

The ionotropic GABAA receptor (GABAAR) has been proven to be an important target of atypical antipsychotics. A novel series of imidazo [1,2-a]-pyridine derivatives, as selective positive allosteric modulators (PAMs) of α1-containing GABAARs with potent antipsychotic activities, have been reported recently. To better clarify the pharmacological essentiality of these PAMs and explore novel antipsychotics hits, three-dimensional quantitative structure–activity relationships (3D-QSAR), molecular docking, pharmacophore modeling, and molecular dynamics (MD) were performed on 33 imidazo [1,2-a]-pyridines. The constructed 3D-QSAR models exhibited good predictive abilities. The dockings results and MD simulations demonstrated that hydrogen bonds, π–π stackings, and hydrophobic interactions play essential roles in the binding of these novel PAMs in the GABAAR binding pocket. Four hit compounds (DS01–04) were then screened out by the combination of the constructed models and computations, including the pharmacophore model, Topomer Search, molecular dockings, ADME/T predictions, and MD simulations. The compounds DS03 and DS04, with higher docking scores and better predicted activities, were also found to be relatively stable in the binding pocket by MD simulations. These results might provide a significant theoretical direction or information for the rational design and development of novel α1-GABAAR PAMs with antipsychotic activities.


Introduction
Schizophrenia, a persistent and chronic psychiatric disorder, affects about 1% of the worldwide population [1,2]. It is mainly characterized by positive symptoms (delusions, auditory hallucinations, and illusions), negative symptoms (anhedonia, apathy, difficulties with concentration, blunted affect, and social dysfunction), and cognitive abnormalities (impairments in working and verbal memory) [3][4][5]. It is reported that the disease ranked third among the most debilitating diseases in the world [6,7]. The treatment of schizophrenia generally focuses on eliminating the disease-associated symptoms. Presently, the main target of most typical antipsychotics is the subcortical dopamine D2 receptor [8][9][10]. Studies have shown that the long-term use of typical antipsychotics might induce side effects, such as oversedation, anesthesia, obesity, myocarditis, liver damage, hyperprolactinemia, postural hypotension, and extrapyramidal symptoms [11,12]. Hence, there is an urgent need to develop novel antipsychotics, engaging different mechanisms of action from those currently used, which would provide alternatives for schizophrenia patients [13,14].
As early as in 1972, GABAergic dysfunction in schizophrenia was first proposed by Roberts [15]. GABA (γ-aminobutyric acid) was first discovered in mammal brains To better explore the crucial pharmacological characteristics of the newly found α1-GABAAR PAMs, and design more-efficient antipsychotic agent leads, herein, 33 imidazo [1,2-a]-pyridine PAMs of GABAARs were selected to perform a systematic modeling study, including three-dimensional quantitative structure-activity relationship (3D-QSAR) models, pharmacophore models, molecular dockings, and molecular dynamics (MD). This study could be clearly illustrated by the following flow chart ( Figure  2). To find novel α1-GABAAR PAMs, virtual screening was then performed based on the best pharmacophore model combined with Topomer Search, molecular dockings, and To better explore the crucial pharmacological characteristics of the newly found α1-GABA A R PAMs, and design more-efficient antipsychotic agent leads, herein, 33 imidazo [1,2-a]-pyridine PAMs of GABA A Rs were selected to perform a systematic modeling study, including three-dimensional quantitative structure-activity relationship (3D-QSAR)

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GA PAMs were selected from the published literatures [26,27]. Their chemical structur activities are given in Table 1. To better understand the QSAR of these novel imidaz a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, inc comparative molecular field analysis (CoMFA), comparative molecular similarity i analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parame all the models are presented in Tables 2 and 3, respectively, using the same traini molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model w employed for the R-group research.

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GABA PAMs were selected from the published literatures [26,27]. Their chemical structures a activities are given in Table 1. To better understand the QSAR of these novel imidazo [1 a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research.

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GABA PAMs were selected from the published literatures [26,27]. Their chemical structures a activities are given in Table 1. To better understand the QSAR of these novel imidazo [1 a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research.

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GABA PAMs were selected from the published literatures [26,27]. Their chemical structures a activities are given in Table 1. To better understand the QSAR of these novel imidazo [1 a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research.

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GABA PAMs were selected from the published literatures [26,27]. Their chemical structures a activities are given in Table 1. To better understand the QSAR of these novel imidazo [1 a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research.

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GABA PAMs were selected from the published literatures [26,27]. Their chemical structures a activities are given in Table 1. To better understand the QSAR of these novel imidazo [1 a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research.

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GABA PAMs were selected from the published literatures [26,27]. Their chemical structures a activities are given in Table 1. To better understand the QSAR of these novel imidazo [1 a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research.

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GABA PAMs were selected from the published literatures [26,27]. Their chemical structures a activities are given in Table 1. To better understand the QSAR of these novel imidazo [1 a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research.

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GABA PAMs were selected from the published literatures [26,27]. Their chemical structures a activities are given in Table 1. To better understand the QSAR of these novel imidazo [1 a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research.

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GABA PAMs were selected from the published literatures [26,27]. Their chemical structures a activities are given in Table 1. To better understand the QSAR of these novel imidazo [1 a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research.

The Statistical Analysis of the 3D-QSAR Models
In the present study, the dataset of 33 imidazo [1,2-a]-pyridines as α1-GABA PAMs were selected from the published literatures [26,27]. Their chemical structures a activities are given in Table 1. To better understand the QSAR of these novel imidazo [1 a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includi comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research. a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includi comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameters all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research. a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameter all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research. a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity indi analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameter all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research. a]-pyridines as α1-GABAAR PAMs, in the present study, the 3D-QSAR models, includ comparative molecular field analysis (CoMFA), comparative molecular similarity ind analysis (CoMSIA), and Topomer CoMFA, were established. The statistical parameter all the models are presented in Tables 2 and 3, respectively, using the same training molecules) and test (8 molecules) sets. Meanwhile, the Topomer CoMFA model was a employed for the R-group research.     Internal and external validation parameters were important criterions for evaluating the quality and credibility of the 3D-QSAR models. In this study, the cross-validation correlation coefficient (q 2 ) value of the CoMFA and Topomer CoMFA models were 0.808 and 0.857, respectively, the non-cross-validated correlation coefficient (R 2 ) value were 0.955 and 0.911, respectively, and the optimum number of components (ONC) were 7 and 15, respectively. These internal validation parameters of the two models satisfied the standards. In the CoMFA model, the contribution rates of the steric (S) and electrostatic (E) fields were 37.3% and 62.7%, respectively, suggesting that the contribution of the E field was more crucial. Additionally, the CoMSIA model could be evaluated by any combination of S, E, hydrophobic (H), hydrogen bond donor (HBD), and hydrogen bond acceptor (HBA) fields. Herein, fourteen different CoMSIA models were built and their statistical results are summarized in Table 2. The CoMSIA-SEHAD model that gave relatively reasonable parameters of internal validation (q 2 = 0.862, R 2 = 0.980, ONC = 13) was chosen as the optimal model for further analysis. In the CoMSIA-SEHAD model, the contributions of the S, E, H, HBA, and HBD fields were 6.8%, 18.0%, 20.8%, 40.7%, and 14.7%, respectively, suggesting that the E, H, and HBD fields played important roles in this CoMSIA model. It was worth mentioning that the r pred 2 values of the CoMFA, Topomer CoMFA, and CoMSIA models were 0.879, 0.935, and 0.927, respectively, indicating that these three models have reasonable predictabilities. All the validation parameters within the standard ranges in Tables 2 and 3 suggested that the established 3D-QSAR models were robust and reliable.
The scatter plots of the actual and predicted activities of all the compounds, by the three models, are shown in Figure 3. The statistical points of all the compounds showed a good linear correlation, which further proved that the 3D-QSAR models have high reliability for predicting the activity of those imidazo [1,2-a]-pyridines as GABA A R PAMs.
contributions of the S, E, H, HBA, and HBD fields were 6.8%, 18.0%, 20.8%, 40.7%, and 14.7%, respectively, suggesting that the E, H, and HBD fields played important roles in this CoMSIA model. It was worth mentioning that the rpred 2 values of the CoMFA, Topomer CoMFA, and CoMSIA models were 0.879, 0.935, and 0.927, respectively, indicating that these three models have reasonable predictabilities. All the validation parameters within the standard ranges in Tables 2 and 3 suggested that the established 3D-QSAR models were robust and reliable.
The scatter plots of the actual and predicted activities of all the compounds, by the three models, are shown in Figure 3. The statistical points of all the compounds showed a good linear correlation, which further proved that the 3D-QSAR models have high reliability for predicting the activity of those imidazo [1,2-a]-pyridines as GABAAR PAMs.

3D-QSAR Contour Map Analysis
The CoMFA, CoMSIA, and Topomer CoMFA contour maps with the most potent compound 14 as a reference are shown in Figure 4-6, respectively. It can be clearly found that the S and E contour maps of the CoMFA models were similar to those of the CoMSIA model. Figure 4, there was a large green contour near the C4 ′ position of the benzene ring, indicating that the introduction of bulky groups at this place might be beneficial for improving the activity. The fact that the methyl or trifluoromethyl substituents at the C4 ′ position were better than the hydrogen atom for the activity could support this result, as illustrated in the following activity orders: 02 (

3D-QSAR Contour Map Analysis
The CoMFA, CoMSIA, and Topomer CoMFA contour maps with the most potent compound 14 as a reference are shown in Figures 4-6, respectively. It can be clearly found that the S and E contour maps of the CoMFA models were similar to those of the CoMSIA model.   As shown in Figure 4, there was a large green contour near the C 4 position of the benzene ring, indicating that the introduction of bulky groups at this place might be beneficial for improving the activity. The fact that the methyl or trifluoromethyl substituents at the C 4 position were better than the hydrogen atom for the activity could support this result, as illustrated in the following activity orders: 02 (CH 3  In the S field contours of the CoMSIA model, a large green contour covering the C 3 position substituent was observed, indicating that bulky volume groups might be more beneficial for the activity at this place. However, there were medium-sized yellow contours around the C 3 position groups, implying that the bulky groups introduced in this area might be disadvantageous for the bioactivity. The volume of the green color block was slightly larger than that of the yellow contour. In consequence, the influence of the green color block on the S field effect in this area was the primary consideration. A similar situation was also presented in the S field of the CoMFA model. Based on the above analysis of the S field color block, we found that further structural modification using bulky substituents at this position might be favorable for the increment of activity. For instance, the introduction of the 3-methylpropanamide group at this position was more beneficial to the activity than the 3-methylacetamide and 3-methylisobutyramide groups. This could be confirmed by the activity order of the following compounds: 14 (propionamide) > 13 (acetamide) > 15 (isobutyramide), 18 (propionamide) > 17 (acetamide) > 19 (isobutyramide), and 22 (propionamide) > 21 (acetamide) > 23 (isobutyramide).
The E field contours of the CoMFA and CoMSIA models are shown in Figure 4. A medium-sized blue contour was directly opposite to the pyrimidine ring, indicating that the introduction of electronegative substituents at this area might be unfavorable to the activity. This might explain why the activity of compound 01 (CH 3 ) was higher than that of compound 07 (CF 3 ). It was also found that the introduction of the F atoms at the C 3 position of the pyrimidine ring had no significant effect on the activity of those PAMs. A red block was close to the C 2 position group, meaning that the electronegative groups at the C 2 position might be advantageous for the activity. This could be verified by the following activity orders: 07 (-F) > 08 (-H) and 10 (-F) > 08 (-H). Meanwhile, a lager red contour was also present near the C 4 position of the benzene ring, indicating that electronegative substituents at this area might be advantageous for the activity. This could be certified by the fact that the activity of compound 10 (-F) was better than that of compound 03 (-H). A medium-sized red contour in the CoMFA model and a lager red contour in the CoMSIA covered up the imidazo [1,2-a]-pyridine ring, indicating that the occupation of electronegative groups at this region was advantageous for increasing the activity. Moreover, a medium-sized red contour near the carboxide on the propionamide group (-CO-), and a blue contour in both the CoMFA and CoMSIA models could be observed lying in the nitrogen atom of the propionamide group (-NH-). This result confirmed that the propanamide group was the essential group to this series of compounds, which was consistent with the S field results.
The H field contours are exhibited in Figure 5. A larger yellow contour appeared on the side of the substituent at the C 4 position of the benzene ring, indicating that hydrophobic groups at this position could be favorable for enhancing the bioactivity. For instance, the hydrogen atom at the C 4 position was replaced by the methyl group, causing the activity of compound 03 (pKi = 6.866) to be less than compound 02 (pKi = 7.602). Meanwhile, in the H field contours of the CoMSIA model, a large yellow contour covering the C 3 position group highlighted the importance of hydrophobic groups in this region. Besides, at the bottom of the C 3 position substituent, the presence of a small white contour indicated that the hydrophobic groups at this region might decrease the bioactivity. We observed that the volume of the yellow color block was significantly larger than the white color block; therefore, we mainly considered the influence of the yellow color block on the hydrophobic effect in this area. The fact that the 3-methylpropionamide substituent at the C 3 position was better for the activity than 3-methylacetamide could explain this result, as illustrated in the following activity orders: 14 (propionamide) > 13 (acetamide), 18 (propionamide) > 17 (acetamide), 22 (propionamide) > 21 (acetamide), 26 (propionamide) > 25 (acetamide), and 30 (propionamide) > 29 (acetamide). The hydrophobicity would become slightly stronger with the length of the alkane chain, which might be compatible with the result of the S field that bulky groups be favorable for enhancing the biological activity at the C 3 position.
The HBD and HBA contours of the CoMSIA model are given in Figure 5. As shown in Figure 6b, a small cyan contour near the C 4 position substituent demonstrated that the presence of HBD groups at this position might favor the activity. This could be corroborated by the following active order: 11 (-F) > 04 (-H). Furthermore, a large cyan contour near the C 3 position substituent of the imidazole ring, and a medium-sized contour were observed near the 3-methylpropionamide substituent on the imidazole ring, demonstrating that HBA groups were favorable for the activity at these positions. These results could well explain why the activities of compounds 14-33 (-CO-) were superior to those of compounds 07-09 (-NH-). In addition, two medium-sized magenta contours near the N atoms of the imidazo [1,2-a]-pyridine ring manifested that the existence of a HBA group on the common scaffold might be significant for the activity. These observations were consistent with a previous study that the 3-methylamide group and the imidazo [1,2-a]-pyridine ring of those PAMs were essential for the bioactivity, as they could serve as HBDs or HBAs to interact with the GABA A R protein [26,27,30].
The contour maps of the S and E fields of Topomer CoMFA ( Figure 6) were basically coincident with that of the CoMFA and CoMSIA models, further verifying the reliability of 3D-QSAR models.
According to the analyses of the 3D-QSAR models, the suitable substituents for increasing the activity of these PAMs at specific regions might be concluded as follows: (1) small group and/or electropositivity at the C 6 position of the imidazo [1,2-a]-pyridine ring; (2) bulky, negative-charged, and/or hydrophobic groups at the C 4 position of the phenyl ring; (3) negatively charged groups at the C 2 position of the phenyl ring; (4) hydrophilic groups at the C 3 position of the imidazo [1,2-a]-pyridine core; (5) 3-methylamide group at the C 3 position of the imidazo [1,2-a]-pyridine core as HBDs or HBAs; (6) HBDs groups at the para-position of the phenyl ring at the C 4 position. Generally, the established 3D-QSAR models were dependable, and could provide a theoretical basis for the design and synthesis of novel efficient PAMs.

Molecular Docking
A newly crystal structure of the human GABA A R α1β2γ2 subtype, in a complex with GABA plus classical benzodiazepines (BZD) diazepam (DZP), has been successfully resolved (PDB ID: 6X3X) in 2020 [32]. In this study, the generated BZD binding pocket of this complex, by the Surflex-Dock method, is depicted in Figure 7a. In order to validate the reliability of the docking method and explore the preferred binding patterns of the test molecules, the co-crystallized DZP was extracted and then re-docked into the BZD binding pocket of the α(+)/γ(−) interface of the human α1β2γ2 GABA A R. As depicted in Figure 7b, the conformation of the re-docked DZP almost overlapped with that of the co-crystallized DZP, and their rotational tendencies were basically the same. The root mean square deviation (RMSD) value of the two conformations was 0.43 Å, which was less than 2.0 Å, illustrating that the binding mode was well-founded and could be applied for further study. The two benzene rings of the DZP molecule were found to form π-π stacking interactions with Tyr210 (loop C of α1 subunit) and Phe77 (loop D of γ2 subunit), respectively. Two benzene rings and the Cl substituent were also involved in hydrophobic interactions with surrounding amino acid residues (Phe100, Tyr160, Val203, Tyr210 in α1, and Tyr58, Phe77 in γ2), which might be one of the main reasons for the stable binding of DZP in the GABA A R protein. These interactions, such as π-π stacking and hydrophobic interactions between DZP and the 6X3X protein, were in accordance with a previous report [33]. The above-mentioned findings further demonstrated that the docking method was accurate and could be used for the following investigation.
To comprehensively survey the binding mechanism of these novel PAMs, 33 imidazo [1,2-a]-pyridine derivatives (Table 1) were then docked into the GABA A R BZD site, using the same method. In general, we found that the docking scores of these PAMs were basically consistent with their activities. The structural models of zolpidem binding at the α(+)/γ(−) interface of 6X3X are displayed in Figure 8a.  Zolpidem was located at the α1/γ2 interface in an oblique insertion, the imidazopyridine ring pointed down in the binding pocket towards His102 (loop A of α1 subunit), and the carbonyl group pointed to loop C. There was a key hydrogen bond formed between the carbonyl oxygen atom and Ser205 in loop C (Ser205-O-H . . . O-C, 1.8 Å). The dimethyl amide side chain embedded into the hydrophobic pocket that was generated by the surrounding amino acid residues, Phe100, His102, Tyr160, Val203, and Tyr210. Additionally, the π-π stacking interaction was also found between the imidazole ring of zolpidem and the benzene ring of the aromatic residues Phe77 (face-to-face) in loop D. All the observed interactions were consistent with the findings reported by Tikhonova et al. [30,34,35]. Compound 14, with the highest potency, and compound 08, with the lowest activity, were selected to concretely analyze the interaction patterns of these PAMs in 6X3X (Figure 8b,c). The docking score of compound 14 (6.91) was higher than that of compound 08 (5.11), which was consistent with their experimental activities. Furthermore, it was noted that compounds zolpidem, 14, and 08 have similar binding orientations and interactions in the binding pocket. Similarly to zolpidem, the imidazole ring, as the common skeleton of compounds 14 and 08, formed a π-π stacking interaction with the aromatic Phe77 (γ1). The dimethylamide side chain of compound 08, and the propanamide moiety of compound 14 were involved in hydrophobic interactions with surrounding amino acid residues. The carbonyl group of compound 14 ( Figure 8b) formed a hydrogen bond with the key residue Ser205 (Ser205-O-H . . . O-C, 1.5 Å). However, this crucial hydrogen bond was not observed in the docking results of compound 08 (Figure 8c). This result was in agreement with the HBA contours of the 3D-QSAR analysis. In addition, the F atom was present at the para-position of the phenyl ring of compound 14, as a HBA formed a hydrogen bond with Ser206 (Ser206-O-H . . . F-Phenyl, 2.6 Å), and this additional weak hydrogen bonding was also presented in the docking result of compound 08. The higher docking score of compound 14 might be due to the stronger hydrogen bond formed with Ser205 (α1), and the deeper binding and rotation of the 3-methylpropionamide part in the hydrophobic pocket.
In comparison with compound 14, compound 08 was observed to lose partial interactions with Ser205 (α1), which might be considered as the main reason that caused its decreased potency. These findings provide evidence to support the essential role of those key residues for the antipsychotic activity in the BZD active site. These results were consistent with the 3D-QSAR analysis, further validating and supplementing the 3D-QSAR models.

Pharmacophore Model
To investigate the key pharmacophore features of these PAMs, in this study, twenty pharmacophore models were generated by aligning and comparing the common features from a test set of ten PAMs, and their statistical results are listed in Table 4. The MODEL_1 was considered to be the optimal pharmacophore model, as it gave the best parameters, including SPECIFICITY = 5.015 (>4), N_HITS = 10, FEATS = 6, PARETO = 0, and EN-ERGY = 25.5. The decoy set method was then used to verify the quality of MODEL_1, and the results are shown in Table S1. These calculated parameters for MODEL_1 could be concluded as follows: GH = 0.734 (0.6 < GH < 1) and EF = 117.398 (>1). These also demonstrated that the MODEL_1 has powerful ability in discriminating the active compounds from the inactive compounds [36]. Therefore, MODEL_1 was selected to analyze the pharmacological features, and was applied for the following virtual screening. Figure 9 displayed six pharmacophore features using the most potent compound 14 as a reference, including three hydrophobic centers (HYs), two hydrogen bond acceptor atoms (AAs), and one hydrogen bond donor atom (DA). The three HYs were located at the center of the imidazole ring, the pyrimidine ring, and the benzene ring, respectively, suggesting the significance of hydrophobic interactions. The common scaffolds of zolpidem analogues were reported to be indispensable to maintain the activity, which was in agreement with the pharmacophore results that three HYs were distributed in the center of the three aromatic rings. Furthermore, we could observe that the HDs were present on the nitrogen atoms of the propionamide moiety, one HA presented on the oxygen atom of the propionamide side chain, and another HA presented on the nitrogen atom of the imidazole ring of the skeleton. These indicated that the propionamide group structure and the imidazole ring might be important for pharmacological action. To sum up, these findings were in accordance with the 3D-QSAR and docking results, further demonstrating the reliability of this pharmacophore model. A graphical SAR summary of these imidazo [1,2-a]-pyridine derivatives, based on the results of the 3D-QSAR models, the molecular dockings, and the optimal pharmacophore model, is shown in Figure 10.

Virtual Screening Analysis
To find potential antipsychotic agent leads, multiple-round virtual screenings were implemented, using the established 3D-QSAR and pharmacophore models in combination with molecular dockings, and the screening workflow is shown in Figure 11. Firstly, the optimal pharmacophore model (MODEL_1) was converted into a UNITY query to screen against the ZINC purchasable database. There were 97,767 compounds that fitted the pharmacophore model, and 4426 compounds with QFIT > 50 were selected for the next round of screening using Topomer Search. In the present work, Topomer Search was employed to screen R groups in the obtained 4426 hit compounds. The screening results were assessed by the TOPDIST and the contribution values of the R groups (TOP-COMFA_R), respectively. Subsequently, 495 R a and 491 R b fragments were gained by Topomer Search, and target fragments with higher contributions were then collected from all the screened R groups using compound 14 as a template. Consequently, 10 R a and 40 R b groups, with the rational TOPDIST and TOPCOMFA_R values, were extracted, and 400 possible compounds were obtained using the rule of permutation and combination. Next, these 400 molecules were docked into the GABA A R BZD site, and 58 compounds (docking score > 9 and Cscore > 4) were obtained. To guarantee that newly designed molecules are the feasible drugs, the ADME/T properties of the 58 newly screened compounds were then predicted by the pkCSM-pharmacokinetics web tool [37], using zolpidem and compound 14 as reference molecules. Their drug-likeness properties were further predicted using the SwissADME online tool [38]. A portion of data was summarized in Table 5. Finally, four newly screened hits (DS01-DS04) were obtained through the following parameter criteria: 150 < MW < 500 g/mol; 20 < TPSA < 130 Å 2 ; high GI absorption; existent blood-brain barrier (BBB) permeability; existent central nervous system (CNS) permeability (log PS > −2); Caco 2 permeability > 0.9; intestinal absorption (human) > 90%; lower inhibitory characteristics with CYP450; a low value of total clearance; synthetic accessibility score < 4; respect all drug-likeness rules; and no skin sensitization.   Table 6 summarizes the chemical structures of four hit compounds (DS01-DS04) with their docking scores and predicted pKi values, by the Topomer CoMFA model. The compounds DS03 and DS04 were selected for further study, as their predicted activities were higher than those of the other hit compounds and compound 14. The docking results of compounds DS03 and DS04 in the BZD active site of human α1β2γ2 GABA A R are depicted in Figure 12. It could be observed that the docking conformations of the compounds DS03 and DS04 were basically similar to that of compound 14.  The important residues within 5 Å of the ligand involved in the interactions are Phe100 (α1), His102 (α1), Tyr160 (α1), Val203 (α1), Ser205 (α1), Ser206 (α1), Tyr210 (α1), Tyr58 (γ2), and Phe77 (γ2). Similarly to zolpidem and compound 14, the carbonyl group in the side chain of compound DS03 formed a hydrogen bond with Ser205 in loop C. Another hydrogen bond was found between the F atom in the para-position of the benzene ring and Ser206. As for compound DS04 (Figure 12b), two hydrogen bonds with residues Ser205 and Ser206 in α1 subunit were also found. Meanwhile, two hit compounds deflected closer to critical residues (e.g., Phe100, His102 in loop A; Tyr160 in loop B; Val203, Tyr210 in loop C) and had multi-hydrophobic interactions with them. In addition, π-π stacking was also generated between the imidazole rings of both compounds and the phenyl groups of the aromatic residue Phe77 (face-to-face) in loop D. A similar π-π stacking interaction was also found between the benzene ring and Tyr58. Therefore, it could be predicted that the hit compounds might have a stronger binding affinity to the 6X3X protein than compound 14. The above-mentioned results revealed that these compounds ought to act as a potential scaffold for designing novel α1-GABA A R PAMs, and might provide meritorious reference for the rational designing of novel antipsychotics leads.

Molecular Dynamics Simulation
In order to further clarify the binding interactions between the screened ligands and the GABA A R protein, the best docking conformations of zolpidem, 14, DS03, and DS04, in the complex with the 6X3X protein, were subjected to 30 ns MD simulations, respectively.
The time-dependent behavior of the complexes 6X3X-DS03 and 6X3X-DS04 in the 30 ns simulation trajectory frames were analyzed using the complexes 6X3X-zolpidem and 6X3X-14 as the comparisons, and their results were summarized in Figure 13. The RMSD analysis usually provides key information about the convergence and stability of these complex systems [37]. We could observe the RMSD values of the protein backbone atoms of four complexes stabilized after 5 ns MD simulations, converging at 0.20 nm (Figure 13a). From the thorough analysis of the trajectory frames, the RMSD values of the backbone atoms of complexes 6X3X-DS03 and 6X3X-DS04 were slightly lower than those of the other complexes. For the ligand RMSD values (Figure 13b), they fluctuated greatly in the 6X3X-zolpidem complex at the beginning of 20 ns, whereas they fluctuated at about 0.06 nm and 0.09 nm in the complexes 6X3X-DS03 and 6X3X-DS04, respectively. Similarly, the RMSD values of four ligands could reach equilibrium finally. Compared to zolpidem and compound 14, compounds DS03 and DS04 had relativity lower RMSD values and fluctuation, indicating a stable interaction or binding to the 6X3X protein. Moreover, the root mean square fluctuation (RMSF) of the Cα residues was computed to identify the structural changes induced by the ligand binding. The RMSF values of the residues in the loops A-F, where the BZD binding site is located, are depicted in Figure 13. From the RMSF trajectories, it is evident that the RMSF values of the most residues in the four complex systems showed similar fluctuations, which demonstrated that four compounds had the analogical binding modes. The RMSF values of the key residues of four complex systems are summarized in Table S2. The critical residues Phe100 (loop A), His102 (loop A), Tyr160 (loop B), Val203 (loop C), Ser205 (loop C), Ser206 (loop C), Tyr210 (loop C), and Phe77 (loop D), in the binding pocket, had relatively low RMSF values (<0.08). These results revealed that the critical residues in the four systems exhibited low flexibility, illustrating that the key interactions of all the compounds in the binding pocket might maintain stability. Moreover, the root mean square fluctuation (RMSF) of the Cα residues was computed to identify the structural changes induced by the ligand binding. The RMSF values of the residues in the loops A-F, where the BZD binding site is located, are depicted in Figure  13. From the RMSF trajectories, it is evident that the RMSF values of the most residues in the four complex systems showed similar fluctuations, which demonstrated that four compounds had the analogical binding modes. The RMSF values of the key residues of four complex systems are summarized in Table S2. The critical residues Phe100 (loop A), His102 (loop A), Tyr160 (loop B), Val203 (loop C), Ser205 (loop C), Ser206 (loop C), Tyr210 (loop C), and Phe77 (loop D), in the binding pocket, had relatively low RMSF values (<0.08). These results revealed that the critical residues in the four systems exhibited low flexibility, illustrating that the key interactions of all the compounds in the binding pocket might maintain stability.  Figure 13e shows that the hydrogen bond numbers between the ligand and protein in the complexes 6X3X-zolpidem, 6X3X-14, 6X3X-DS03, and 6X3X-DS04, during the MD simulation, fluctuated 0-1, 0-3, 0-4, and 0-4, respectively. It is worth noting that the hydrogen bond numbers of the complexes 6X3X-DS03 and 6X3X-DS04 were slightly higher than that of the complexes 6X3X-14 and 6X3X-zolpidem. Furthermore, the hydrogen bond numbers in the complexes 6X3X-zolpidem, 6X3X-14, 6X3X-DS03, and 6X3X-DS04 maintained in 1, 1, 2, and 3 after 10 ns, respectively. This indicated that the hit compounds DS03 and DS04 might be more stable than zolpidem and compound 14 in the binding pocket. Afterwards, we estimated the gyration radius (Rg) of the Cα atoms of the four complex systems, to evaluate the compactness of the protein structure during MD simulations (Figure 13f). All the complexes showed slight fluctuations in the first 20 ns, and the Rg values remained around 3.85 nm finally, suggesting that the protein conformations might be basically stable in the 30 ns dynamic simulation, which was consistent with the RMSD results.
The molecular mechanics-Poisson Boltzmann surface area (MM-PBSA) approach was considered to be one of the most suitable procedures to calculate the binding free energies (∆G binding ), which has high accuracy and high computational potency in calculating the binding affinities of ligands with their targets [39]. As shown in Table 7, the ∆G binding values of the compounds zolpidem, 14, DS03, and DS04 in the GABA A R protein were −95.181, −120.055, −131.507, and −126.376 kJ/mol, respectively, which was consistent with the docking results. The van der Waals energy was observed to be the largest contributor to ligand binding, whereas the polar solvation energy was disadvantageous for the binding.
In conclusion, zolpidem, compound 14, DS03, and DS04 could be stable in binding to the 6X3X protein during the whole MD simulation, and the screened hits DS03 and DS04 might have better binding stability than zolpidem and compound 14. Table 7. The free binding energies of zolpidem, compound 14, and the newly screened hit compounds (DS03 and DS04) in the human GABA A R (PDB ID: 6X3X).
( Figure 14a). Topomer CoMFA was a fragment-based 3D-QSAR method that integrates the initial "topomer" methodology and the CoMFA technology [43]. Unlike CoMFA and CoMSIA, the most critical step was to select the split mode or identify the R group in Topomer CoMFA modeling. The prediction accuracy of the Topomer CoMFA model and the reliability of the contour map depend on this step [44]. Each molecule can be divided into several small fragments by cutting the twistable chemical bond, and then the S and E fields of the corresponding fragments could be calculated, respectively. In this study, the Topomer CoMFA model was constructed by using the same training and test sets of the CoMFA model. The Topomer CoMFA model was constructed by dividing compound 14 into two segments, R a (magenta) and R b (blue), as shown in Figure 14. The segmentation position was shown as a black curve in Figure 14b, and the other training set compounds were automatically identified and cut in the same pattern. Subsequently, three-dimensional conformations of the segments were obtained and used for predicting their biological activities or properties.

Analysis and Validation of the QSAR Model
In order to evaluate the reliability and predictive ability of the established 3D-QSAR model, internal and external verification parameters were subsequently statistically analyzed [36]. The q 2 value and the ONC were obtained by means of leave-one-out crossvalidation analysis. The R 2 , the Fisher's statistic values (F), and the standard error of estimate (SEE), were calculated to assess the predictive ability of models. The abovementioned parameters were usually considered as internal validation parameters. A model, which is equipped with the following internal parameter ranges: q 2 > 0.5, R 2 > 0.6, and SEE << 1, could be considered as a trustworthy model and might have good internal prediction capabilities [45]. Furthermore, external verification parameters, including r 0 2 , r 0 2 , k, k , r m 2 , r m 2 , r pred 2 , ∆r m 2 , r m 2 and RMSE, were further taken into consideration [41]. The QSAR model would be deemed to have excellent external predictive ability only if those statistical parameters met following requirements: r pred 2 > 0.6, r 0 2 (r 2 − r 0 2 )/r 2 < 0.1 or r 0 2 (r 2 − r 0 2 )/r 2 < 0.1, 0.85 < k < 1.15 or 0.85 < k < 1.15, ∆r m 2 < 0.2, r m 2 , and > 0.5 [36,45].

Molecular Docking
In this study, the molecular docking was performed using the SYBYL-X 2.1 software and was visualized using PyMol 2.3.3 software (DeLano Scientific LLC, San Carlos, CA, USA). The newly resolved crystal structure of human α1β2γ2 GABA A R (PDB ID: 6X3X) with a resolution of 2.92 Å was used as a receptor, which contained a co-crystallized ligand DZP at the BZD binding site [32]. After the pretreatment steps of the original protein, including hydrogenating, adding electric charges, fixing protein side chains and termini chains, removing waters, and extracting the co-crystallized ligand, the applicable binding pocket was generated by a ligand mode. In order to examine the dependability of the docking method, the extracted DZP was firstly redocked into the generated binding pocket by the Surflex-Dock Geom module to inspect whether the redocked ligand conformation and the originally crystallographic conformation overlap. Meanwhile, the conformation differences between the redocked and original ligands were estimated by the RMSD value. An RMSD < 2.0 Å was regarded as a reference standard, suggesting that the used docking method is credible [41]. After that, all selected imidazo [1,2-a]-pyridine PAMs were docked into the active pocket using the same method.

Pharmacophore Model
A pharmacophore is a molecular interaction characteristic shared by a group of active molecules. In this study, ten imidazo [1,2-a]-pyridine as GABA A R PAMs (Table 1) with relatively high biological activities and diverse structures were selected to establish pharmacophore models by the genetic algorithm with linear assignment of hyper-molecular alignment of datasets (GALAHAD) module of SYBYL-X 2.1 [46]. The remaining compounds were used for the evaluation of the constructed models. Twenty models were generated, and the model with high values of SPECIFICITY, N_HITS, HBOND, MOL_QRY and the low value of ENERGY was selected for further study. To confirm whether the chosen model was sufficient to recognize the active compounds from the database, it was necessary to evaluate the model using the decoy set method [47]. The potential pharmacophore models were performed to retrieve a decoy set database, which was composed of 3892 inactive compounds (downloaded from http://dud.docking.org/r2/, accessed on 1 October 2020) and 23 active compounds (Table 1) in addition to those used to generate this pharmacophore model. The following two main indexes were used to assess the model: the enrichment factor (EF) and the Güner-Henry (GH) score, defined by Equations (1) and (2), respectively.
Ha, Ht, A, and D represent the number of true positive compounds in the database, the number of all hit compounds, the number of true positive compounds in the list, and the total number of compounds in the decoy database, respectively. Generally, a reliable model needs to satisfy EF > 1 and 0.6 < GH < 1 [36,47].

Virtual Screening
Considering structural diversity and commercial availability, the ZINC purchase database (http://zinc15.docking.org, accessed on 1 November 2020) that contained about 20 million compounds was used for virtual screening in this study. A multi-stage virtual screening was carried out against the database through the combination of the optimal pharmacophore model, the Topomer Search technology, molecular dockings, and ADME/T predictions. In the first stage, the extraction of pharmacophore features from the best pharmacophore model was used as a 3D search query to retrieve potential molecules. The QFIT parameters were used to evaluate the matching degree of screened compounds with the pharmacophore features. Then, the compounds with a QFIT > 50 were selected for the second-round screening of Topomer search. The Topomer Search technology decomposes all database molecules into different groups, and compares topological similarities with the suspicious R groups obtained from the Topomer CoMFA model. Subsequently, the contribution of a series of R groups to the activity was predicted by the established Topomer CoMFA model. In this study, the R a and R b segments were selected as templates ( Figure 14). The corresponding R groups were screened from the database, and the TOPDIST and the minimum heavy atom of each fragment were set to 185 and 3, respectively, to evaluate the degree of binding [36]. Target fragments with higher contributions were collected from all screened R groups using compound 14 as a template.
Following, some newly molecules were designed from the combination of the hit R a and R b groups for the next round of docking screening. The hit compounds with a docking score of >9 and good binding patterns were selected for further study. Finally, ADME/T properties of identified hit compounds, compound 14, and zolpidem were predicted by two web tools of pkCSM-pharmacokinetics (http://biosig.unimelb.edu.au/pkcsm/prediction, accessed on 1 November 2020) [37] and SwissADME (http://www.swissadme.cn, accessed on 1 November 2020) [38]. Following properties were considered preferentially to obtain the satisfactory compounds: physicochemical properties, lipophilicity, water solubility, blood-brain barrier permeability (logBB), CNS permeability, synthetic accessibility, and toxicity. The hit compounds with desired pharmacophore features and Topomer CoMFA compliance, high docking scores, good binding modes, and ideal ADME/T evaluation results were further studied by MD simulations.

Molecular Dynamics
To further confirm the stability of the hit compounds in the human GABA A R, 30 ns MD simulations were performed on different protein complexes using GROMACS 2016.5 software (Uppsala University, Stockholm University, and the Royal Institute of Technology, Sweden). The protein topological parameters were created by the Pdb2gmx module under the AMBER99SB force field [48]. The ligand topological files were automatically generated by the ACPYPE tools. The molecular system was solvated with SPC/E water models with about a 12 Å buffering distance between the protein receptor and the edges of the cubic box, and twenty additional chloride ions for the neutralization of the system. Moreover, the protein-ligand complexes were successively subjected to 100 ps simulations to accomplish the NVT and NPT equilibrium at 300 K and 1 atm [49]. Finally, a 30 ns MD simulation by 2 fs per step was carried out for each complex. A few parameters, such as RMSD, RMSF, Rg and hydrogen bond numbers, were conducted to investigate the stability or variation in all complex systems during the 30 ns dynamic simulation. The binding free energy was calculated by using the MM-PBSA method [39], and the calculation formula of binding free energy was calculated as follows: ∆G bind = G complex − G f ree−protein − G f ree−ligand (3) in which G complex represents the complex energy, G free-protein represents the receptor energy and G free-ligand represents the energy of the unbound ligand. Binding energies were extracted using the MmPbSaStat.py program and key residue contributions towards binding were gained.

Conclusions
In this study, an integrated set of computational methodologies was used to explore the novel hit compounds of α1-GABA A R PAMs. The 3D-QSAR models were constructed based on 33 imidazo [1,2-a]-pyridines, to visually understand the effect of diverse substitutions on the activity of these GABA A R PAMs. The molecular docking studies revealed the interaction patterns of these PAMs in the BZD pocket, indicating that these imidazo [1,2-a]pyridines could form key hydrogen bonds with the residues Ser205 (α1) and Ser206 (α1), and could form hydrophobic or π-π stacking interactions with the residues Phe100 (α1), His102 (α1), Tyr160 (α1), Val203 (α1), Tyr210 (α1), and Phe77 (γ2). These interactions might be essential for the bindings or activities of these GABA A R PAMs. The pharmacophore model further confirmed that the hydrophobic aromatic ring and the side chain of the acylamino group were important pharmacophore features. The best pharmacophore model, containing six key features, was in accordance with the 3D-QSAR models and docking results. The virtual screening based on the best pharmacophore model, Topomer Search, molecular dockings, ADME/T predictions, and MD simulations finally confirmed four newly hit compounds (DS01-DS04). The compounds DS03 and DS04 could steadily bind to the 6X3X protein in the dynamic environment, and this was further corroborated by the binding energy analysis; however, its synthesis and biological activity remain to be further studied. We believe that these findings could provide profound theoretical direction and information for the further design and development of novel α1-GABA A R PAMs with good antipsychotic activities.