Use of connectivity index and simple topological parameters for estimating the inhibition potency of acetylcholinesterase.

Acetylcholinesterase (AChE) has proven to be an effective drug target in the treatment of neurodegenerative diseases such as Alzheimer's, Parkinson's and dementia. We developed a novel QSAR regression model for estimating potency to inhibit AChE, pK i, on a set of 75 structurally different compounds including oximes, N-hydroxyiminoacetamides, 4-aminoquinolines and flavonoids. Although the model included only three simple descriptors, the valence molecular connectivity index of the zero-order, 0 χv , the number of 10-membered rings (nR10) and the number of hydroxyl groups (nOH), it yielded excellent statistics (r = 0.937, S.E. = 0.51). The stability of the model was evaluated when an initial set of 75 compounds was broadened to 165 compounds in total, with the increase of the range of pK i (exp) from 6.0 to 10.2, yielding r = 0.882 and S.E. = 0.89. The predictive power of the model was evaluated by calculating pK i values for 55 randomly chosen compounds (S.E.test = 0.90) from the calibration model created on other 110 compounds (S.E. = 0.89), all taken from the pool of 165 compounds.

Introduction

Acetylcholinesterase (AChE) have proven to be effective in the treatment of Alzheimer’s and Parkinson’s disease symptoms. The current treatment is based on AChE inhibitors including donepezil, rivastigmine and alkaloid galantamine (Giacobini, 2006, Mohammad et al., 2017, Xie et al., 2020). Although Tacrine (Cognexw®) was approved as a drug for AD treatment, it was discontinued from medical use due to high hepatotoxicity. Therefore, with the ageing of the world population and increased risk of dementia, the development of AChE inhibitors attracts the highest scientific interest in the process of designing safer and more effective drugs (Sanad and Mekky, 2021, Xie et al., 2020).

The QSAR (quantitative structure–activity relationship) method represents an important tool for drug development and has led to numerous AChE QSAR models (Jana et al., 2018, Kumar et al., 2020, Niu et al., 2017) of different complexity and predictivity. The development of a QSAR regression model could facilitate the development of therapeutic ligands by establishing a correlation between the chemical functionalities of the ligand and the desired biological activity. The proposed regression model comprised of the molecular parameters (descriptors) of interest would enable the prediction of biological activity and ease the design of new compounds with the desired activity (Kubinyi, 1993, Karelson, 2000, Selassie and Verma, 2010). The activity prediction of a QSAR model and its accuracy is based on the selection of appropriate molecular descriptors and the reliability of the measured biological activity (Leach, 1996, Shityakov et al., 2014).

The potential of QSAR models using scoring functions to predict the inhibition potency of acetylcholinesterase (AChE) ligands was analyzed in a recent study (Šinko, 2019). The study indicated that the PLP2 scoring function predicts the inhibition potency of ligands with a coefficient of determination r2 = 0.591. Several scoring functions were tested against AChE-ligand complexes deposited in the PDB base: LigScore1, LigScore2, PLP1, PLP2, Jain, PMF and PMF04. The study showed that the drawback of the scoring function evaluation was the low uniformity of kinetic data (Ki or IC50) obtained using various methods of determination and the enzyme source. Kinetic data were collected under different experimental conditions, i.e. temperature, as well as using various enzyme species as a source of AChE, one species for data measurements and another for the determination of the crystal structure of the AChE-ligand complex.

It was interesting to see that studies using the same ligand, e.g. galantamine, and AChE from the same source Electrophorus electricus may differ in the obtained results due to measuring the IC50 (0.36–1.07 μM) (Atanasova et al., 2015, Mary et al., 1998) instead of the Ki (0.19 μM) (Rahman et al., 2006). Uncertainty in the IC50 value determination is caused by the type and concentration of the substrate used for measurements and therefore Ki is a more reliable parameter, as it is a measure of enzyme ligand affinity in the absence of a substrate. The negative effect of the solvent mixture on the AChE enzyme activity when using ethanol or DMSO buffer should be tested, as ethanol or DMSO apparently increase the ligand inhibition potency due to AChE inhibition (Fekonja et al., 2007, Kumar and Darreh-Shori, 2017). Therefore, the effect of a solvent on AChE activity needs to be characterized and compensated properly.

In a study by Wong et al., the problem of various enzyme sources in a QSAR analysis of tacrine-like inhibitors was reported. Therefore, they created 10 different QSAR regression models for each AChE source, e.g. human, Electrophorus electricus and bovine AChE (Wong et al., 2014). The inhibition potency of tacrine-like inhibitors was evaluated by the Ellman (Ellman et al., 1961) or Rappaport method (Rappaport et al., 1959) using acetylthiocholine iodide or acetylcholine chloride, respectively, as the substrate. The reported IC50 values, obtained using different experimental conditions, led to the development of different QSAR regression models to increase model predictivity and overcome the problem of experimental conditions.

The goal of this study was to highlight the key structural features of AChE ligands in correlation with ligand pKi values indicating inhibitory activity using the simplest possible QSAR model. The different QSAR models presented in the literature describe various parameters for AChE ligands, but often these parameters cannot be easily linked to a ligand’s physicochemical properties (Gurung et al., 2017, Jana et al., 2018, Šinko, 2019, Wong et al., 2014). In our approach, we developed a simple QSAR model for the prediction of the human AChE inhibition constant, pKi, on a set of 75 compounds including 4-aminoquinolines, oximes, flavonoids and N-hydroxyiminoacetamides. For all of the compounds in this study, Ki was measured by our laboratory and published previously (Tables 1 and S3, Fig. 1) (Bosak et al., 2019, Bosak et al., 2017, Bušić et al., 2016, Katalinić et al., 2010, Kovarik et al., 2008, Maček Hrvat et al., 2020, Maraković et al., 2020, Maraković et al., 2016, Šinko et al., 2010, Zandona et al., 2020). The assay used for the AChE activity measurement was based on the Ellman method (Ellman et al., 1961), with standardized activity measurement regarding enzyme, substrate and inhibitor concentrations (Eyer et al., 2003, Reiner et al., 2000). To avoid the artefacts of Ki calculation, AChE inhibition was limited for 20–80% of the control activity (Bosak et al., 2019, Simeon-Rudolf et al., 2001). Moreover, AChE activities were corrected when the oxime-induced degradation of the substrate (oximolysis) was above 10% of the enzyme control activity (Maček Hrvat et al., 2018, Šinko et al., 2007, Šinko et al., 2006).

Table 1

The values of negative logarithms of the AChE inhibition constant (pKi), and molecular descriptors for 75 compounds. 0χv, nOH and nR10 were calculated by the E-DRAGON program system. The compound names are the same as in the original papers whose references are given.

No.	Compound	pKi	0χv	nOH	nR10
1	ICD-585a	4.55	11.70	1	0
2	HI-6a	4.51	11.40	1	0
3	HLo-7a	4.62	12.79	2	0
4	K027a	4.14	11.70	1	0
5	K048a	3.96	12.41	1	0
6	K033a	4.77	12.39	2	0
7	TMB-4a	3.74	11.68	2	0
8	DMB-4a	4.00	10.98	2	0
9	MMB-4a	3.31	10.27	2	0
10	ICD-692a	4.74	11.99	1	0
11	ICD-467a	5.92	11.53	1	0
12	K114b	5.68	14.29	2	0
13	K127b	3.76	12.81	1	0
14	K203b	4.44	12.15	1	0
15	Ic	2.93	9.610	1	0
16	IIc	3.45	14.08	1	0
17	IIIc	4.31	14.93	1	0
18	IVc	3.87	15.41	2	0
19	CQd	4.96	7.050	0	1
20	CQ2d	5.39	8.96	0	1
21	TFCQ2d	5.44	9.46	0	1
22	TFCQ8d	6.34	13.70	0	1
23	CQ8d	6.21	13.20	0	1
24	CQAdd	6.11	16.77	0	1
25	Chloroquined	5.40	14.53	0	1
26	CQEtOHd	5.00	8.83	1	1
27	Galangine	4.07	10.20	3	1
28	Kaempferole	4.03	10.57	4	1
29	Quercitine	4.42	10.94	5	1
30	Myricetine	4.42	11.31	6	1
31	Luteoline	4.18	10.57	4	1
32	Fisetine	4.00	10.57	4	1
33	Apigenine	3.92	10.20	3	1
34	Rutine	3.52	22.29	10	1
35	Metaproterenolf	2.51	8.94	3	0
36	Terbutalinef	2.32	9.86	3	0
37	Fenoterolf	3.07	12.40	4	0
38	Epinephrinef	2.19	7.36	3	0
39	Isoproterenolf	2.60	8.94	3	0
40	Isoetharinef	3.68	10.51	3	0
41	Salbutamolf	2.70	10.57	3	0
42	Salmeterolf	4.52	17.93	3	0
43	1g	3.64	12.15	3	0
44	2g	3.36	12.45	3	0
45	3g	3.90	13.21	3	0
46	4g	4.12	14.04	3	0
47	5g	3.78	13.34	3	0
48	6g	3.98	13.08	3	0
49	7g	4.05	13.48	3	0
50	8g	3.90	13.48	3	0
51	9g	3.90	15.46	3	0
52	1ah	7.82	23.95	1	1
53	2ah	8.22	24.65	1	1
54	1bh	8.05	24.65	1	1
55	2bh	7.55	25.36	1	1
56	1ch	7.17	25.36	1	1
57	2ch	7.49	26.07	1	1
58	1dh	6.89	23.95	1	1
59	2dh	7.64	24.65	1	1
60	1eh	7.39	24.65	1	1
61	2eh	7.00	25.36	1	1
62	Q1i	2.42	5.96	1	0
63	Q2i	3.12	7.01	1	0
64	Q3i	3.28	8.00	1	0
65	Q4i	3.26	8.71	1	0
66	Q5i	5.22	14.79	1	0
67	Q6i	3.80	10.1	1	0
68	Q7i	4.49	11.99	1	0
69	Q8i	4.17	11.99	1	0
70	Q9i	4.39	11.16	1	0
71	Q10i	4.04	11.16	1	0
72	Q11i	4.57	11.02	1	0
73	Q12i	3.87	11.02	1	0
74	Q13i	4.74	11.29	1	0
75	Q14i	4.31	11.29	1	0

From Ref. (Kovarik et al., 2008).

From Ref. (Šinko et al., 2010).

From Ref. (Maraković et al., 2016).

From Ref. (Bosak et al., 2019).

From Ref. (Katalinić et al., 2010).

From Ref. (Bosak et al., 2017).

From Ref. (Bušić et al., 2016).

From Ref. (Maček Hrvat et al., 2020).

From Ref. (Zandona et al., 2020).

Fig. 1

Example of chemical structures of the 10 groups of compounds used in the set of 75 compounds.

Materials and methods

Calculation of topological indices

Molecular descriptors were calculated by the E-DRAGON program developed by Tetko et al. E-DRAGON provides more than 1 600 molecular descriptors (topological, constitutional, geometrical, etc.) in a single run (Tetko et al., 2005). The connectivity matrices were constructed using the Online SMILES Translator and Structure File Generator (Online SMILES Translator and Structure File Generator, 2020). The SMILE formulas, for all compounds studied, are given in Supplement (Tables S3, S4 and S5).

The model developed in this study is based on the topological 0χ index (the valence molecular connectivity index of the zero-order) (Kier and Hall, 1986, Kier and Hall, 1976a, Kier and Hall, 1976b, Randić, 2008), which was defined as:where δ(i) is the weight (valence value) of each vertex (atom) i in a vertex-weighted molecular graph. The valence value, δ(i), of vertex i is defined as:where Z(i) is the number of valence electrons belonging to the atom corresponding to vertex i, Z(i) is its atomic number, and H(i) is the number of hydrogen atoms attached to it. For instance, the delta values for the primary, secondary, tertiary, and quaternary carbon atoms are 1, 2, 3, and 4, respectively, while for the oxygen in the OH group, this equals 5 and for the NH2 group 3. It should be pointed out that 0χ is the only one of the many members from the family of valence connectivity indices , which differ amongst each other by path length, i.e. the number of consecutive chemical bonds. From Eq. (1) it can be seen that 0χ has a path order of zero, i.e. it considers only separate vertices (atoms). 1χ (1χ = Σ [δ(i) δ(j)]−0.5) considers vertices (atoms) i and j, making up a path with a length of 1 (one consecutive chemical bond), 2χ (2χ = Σ [δ(i) δ(j) δ(k)]−0.5) considers vertices (atoms) i, j and k, making up a path with a length of 2 (two consecutive chemical bonds), etc. Connectivity indices are also called branching indices and are among the most used topological indices in QSPR/QSAR, e.g. 3χ was very successfully used for the estimation of the stability constants of metal chelates (Miličević and Raos, 2008; Raos et al., 2008).

Regression calculations

Regression calculations, including the leave-one-out procedure (LOO) of cross validation, were done using the CROMRsel program (Lučić and Trinajstić, 1999). The standard error of the cross-validation estimate was defined as:where ΔX and N denote cv residuals and the number of reference points, respectively.

Results and discussion

Although the correlation of the valence molecular connectivity index of the zero-order, 0χ on pKi, yields somewhat worse statistics (r = 0.795 and S.E. = 0.88, N = 75) than the correlation with the squared Ghose-Crippen octanol–water partition coefficient, AlogP2, (r = 0.857 and S.E. = 0.74, N = 75) and a few other topological indices, it has captured our attention. More precisely, considering the presence of a 10-membered ring (two fused six-membered rings) in the molecule, two almost parallel correlation lines on 0χ pKi dependence can be drawn (Fig. 2). The first line (triangles in Fig. 2, r = 0.800, N = 26) belongs to molecules with a 10-membered ring and the second (circles in Fig. 2, r = 0.575, N = 49) is without a 10-membered ring. It can also be seen that molecules with a 10-membered ring in their structure generally have higher values of pKi. Flavonoids (27–34) are the only compounds for which this does not apply, but their structure is highly rigid in comparison with other compounds with a 10-membered ring. This is especially true for the flavonoid rutin (34), which by far has the highest number of OH groups (10 hydroxyl groups) of all molecules in the set.

Fig. 2

Dependence of pKi on the 0χ index for 75 molecules in the set (Table 1). The lines of correlations are made on two subsets of molecules, with (N = 26, r = 0.800) and without (N = 49, r = 0.575) a 10-membered ring. Triangles denote molecules with a 10-membered ring.

All these were the reason why we added nR10 and nOH descriptors alongside 0χ into the equation. In that way, we developed three descriptor model for the estimation of pKi:yielding r = 0.937, S.E. = 0.51 and S.E.cv = 0.53 (a = 2.54(19), b1 = 0.170(13), b2 = 1.13(14), b3 = –0.353(38)) for the set of 75 compounds, Fig. 3. It is also important to note that the correlations between the pairs of descriptors were very small; r = 0.233, 0.123 and 0.483 for nOH vs. nR10, 0χ nOH and 0χ nR10, respectively.

Fig. 3

Plot of experimental vs. calculated (using Eq. (4)) pKi values; N = 75, r = 0.937, S.E. = 0.51 and S.E.cv = 0.53.

Some topological and constitutional descriptors correlated to pKi showed similar statistics and a similar pattern as 0χ, like the valence molecular connectivity index of the first order, 1χ, the eccentric connectivity index, CSI, and the number of atoms, nAT (r = 0.809, 0.798 and 0.777, respectively). Their implementation in Eq. (4) in place of 0χ yielded slightly worse statistics (S.E. = 0.53, 0.55 and 0.54, respectively) than the standard error obtained by 0χ (S.E. = 0.51).

Although the best possible model with three descriptors chosen among all of the 1399 calculated descriptors gave better results than Eq. (4) (r = 0.952, S.E. = 0.44 and S.E.cv = 0.47), the descriptors used in that model were not easy to connect to the structure of compounds; highest eigenvalue number of Burden matrix weighted by atomic Sanderson electronegativities (BEHe1), 3D-MoRSE - signal 13 weighted by atomic van der Waals volumes (Mor13v) and the difference between multiple path count and path count (PCD).

Previously (Šinko, 2019) we evaluated models using scoring functions for the pKi (or pIC50) estimation of 56 molecules (Tables S1 and S4). By applying our model (Eq. (4)) on the same set of compounds, the statistics were not so good, r = 0.830, rcv = 0.798, S.E. = 1.20 and S.E.cv = 1.30, but one must be aware that the Ki (AChE) values for this set were not measured by the same laboratory and on the same type of AChE (they used human, mouse, etc.). Moreover, instead of Ki, for some molecules IC50 values (Rahman et al., 2006, Atanasova et al., 2015, Herkert et al., 2011, Mary et al., 1998, Saxena et al., 1999) were given. However, when we brought together this set of 56 compounds with our set of 75 compounds, the results of regression on 131 molecules were very good (N = 131, r = 0.892, rcv = 0.883, S.E. = 0.94 and S.E.cv = 0.97), especially as the range of experimental pKi (or pIC50) increased from 6.03 to 10.21.

We also used 34 oximes, Tables S2 and S5, from our previous paper (Katalinić et al., 2016), where we showed that pIC50 can be estimated by using only one topological index; the model using eccentric connectivity index (CSI, Sharma et al., 1997) yielded r = 0.957, S.E. = 0.21 and S.E.cv = 0.23. The 0χ index yielded slightly worse but also excellent results (r = 0.926, S.E. = 0.27 and S.E.cv = 0.29). Combining this set of 34 compounds with two sets presented above (N = 75 and N = 56), our model (Eq. (4)) yielded r = 0.882, S.E. = 0.89 and S.E.cv = 0.91 (a = 3.17(23), b1 = 0.172(17), b2 = 0.753(90), b3 = –0.451(55); N = 165, Fig. 4). Comparing these regression parameters (on a set of 165 molecules) with the parameters in the model for the 75-member set (a = 2.54(19), b1 = 0.170(13), b2 = 1.13(14), b3 = –0.353(38)) great similarity can be observed, although within the limits of S.E. only for 0χ (b1 = 0.172(16) vs. 0.170(13)).

Fig. 4

Plot of experimental vs. calculated (using Eq. (4)) pKi (or pIC50) values; N = 165, r = 0.882, rcv = 0.874, S.E. = 0.89 and S.E.cv = 0.91. Circles denote the set of 75 compounds, triangles the set of 56 compounds used in our previous report (Šinko, 2019), and empty circles the set of 34 oximes (Katalinić et al., 2016).

We tested the predictability of our model (Eq. (4)) by a training/test method. We selected every third molecule (molecules 3, 6, 9, 12…, Table 1, S1 and S2) into the test set and thereby divided the set of 165 compounds into a training set (110 compounds) and test set (55 compounds). Statistics of the calibration model calculated from Eq. (4) on the training set (r = 0.882, S.E. = 0.89 and S.E.cv = 0.92, N = 110) were of the same quality as the model made on 165 compounds and we used it for predicting the pKi values of 55 molecules from the test set. The standard error of the test set (S.E.test = 0.90) was very similar to the S.E. and S.E.cv of the calibration model (S.E. = 0.89 and S.E.cv = 0.92), which proved the high predictive power of Eq. (4).

A comparison of the AChE active site amino acid composition and related functional characteristics with the QSAR descriptors 0χ, nR10 and nOH led us to the following observations. The human AChE active site gorge is a ∼20 Å deep and ∼5 Å wide cavity composed of mainly aromatic residues (Phe, Trp or Tyr) thus creating a hydrophobic space (Ordentlich et al., 1993, Sussman et al., 1991). At the bottom of the narrow active site, where substrate hydrolysis occurs, a catalytic triad Ser203, Glu334 and His447 is located (Fig. 5). The substrate of AChE is a small carboxyl ester with a positively charged choline part, acetylcholine. During acetylcholine hydrolysis, the following key interactions between enzyme residues and substrate formed: hydrogen bonds, hydrophobic interactions and cation-π interactions (Colletier et al., 2006). Ligands that can create these interactions producing strong binding within the AChE active site are possible drug candidates.

Fig. 5

Crystal structure of the human AChE active site (PDB ID 4EY7). Aromatic residues of the peripheral anionic site are orange, of the choline binding site purple, and of the acyl pocket blue. The anti-Alzheimer drug donepezil is yellow. Catalytic Ser203 is represented by a ball and stick. Figure adopted from ref. (Miličević and Šinko, 2021).

Several residues of the AChE active site: Asp74, Glu202, Tyr124, Ser293 and Tyr337 have hydrogen bond donor or acceptor groups, and therefore may stabilize ligands via hydrogen bonds (Šinko, 2019). Hydrogen bond donors or acceptor groups are the molecular basis for an nOH descriptor presence in the QSAR model. Two important sub-domains of the AChE catalytic site, the peripheral anionic site and choline binding site, are responsible for substrate transport and orientation during catalytic turnover (Colletier et al., 2006). Tryptophan Trp86 and Trp286 are key residues of the choline binding site and peripheral anionic site, respectively (Colletier et al., 2006, Ordentlich et al., 1993). These two residues create cation-π and/or π-π interactions with ligands having aromatic groups; e.g. two fused benzene rings (nR10) in quinolines or benzopyrans create π-π stabilizing interactions by overlapping with Trp indole ring. Studies have shown that ligands long enough to bind simultaneously in the choline binding site and the peripheral anionic site are more potent inhibitors of AChE, e.g. donepezil depicted in Fig. 5 (Bourne et al., 2016, Cheung et al., 2012, Felder et al., 2002, Rydberg et al., 2006). Crystal structures of AChE-ligand complexes showed that the binding would be stronger if the conformation of the residues stabilizing the ligand was complementary to the conformation of the apo-AChE active site, e.g. donepezil only affects the conformation of Tyr337, while conformation of the other residues remained unchanged (Cheung et al., 2012, Gerlits et al., 2019). It has been shown that some ligands upon binding induce the change of the conformation of the AChE active site residues (Bourne et al., 2010), therefore there is no clear connection between active site residues conformation and the strength of ligand binding (Šinko, 2019). Index 0χ accounts for the size of the ligand but more importantly for the complexity of the ligand structure, including branching. We showed a positive correlation between 0χ and pKi in Fig. 2. A similar finding was presented in a prior study where the complexity property of ligands positively correlated with the scoring functions was described (Šinko, 2019).

Conclusions

The presented QSAR model (Eq. (4)) is based on the valence molecular connectivity index of the zero-order, 0χ, combined with the number of 10-membered rings (nR10) and the total number of OH groups in a molecule (nOH). On a set of 75 molecules, the model yielded S.E. = 0.51, meaning that pKi (or pIC50) values can be estimated by an error of 8.5% of the pKi range ((S.E./range pKi) 100%). Although the model on 131 compounds, after adding 56 compounds from the literature (Rahman et al., 2006, Atanasova et al., 2015, Herkert et al., 2011, Mary et al., 1998, Saxena et al., 1999) seemed to deteriorate (S.E. = 0.94), ultimately this was not the case. The reason was that the pKi (or pIC50) range increased from 6.03 to 10.21, so, although the S.E. almost doubled, the error of estimation increased only slightly, to 9.2%. This, and especially the regression on 56 molecules (S.E. = 1.20, with an error of estimation of 12.8%), told us that QSAR should be avoided on non-standardized experimental data. When we added 34 oximes (IC50 measured in our laboratory) to the set of 131 compounds, the range of pKi (or pIC50) values stayed the same, and the S.E. and error of estimation dropped to 0.89 and 8.7%, respectively.

Comparing errors of estimation yielded by Eq. (4) with the mean experimental error in Ki measurements, which was 15% for the set of 75 molecules, we can conclude that our results are very satisfactory (Raos et al., 2008, Raos and Miličević, 2016). This is proof that the variables we used in our three-descriptor model (0χ, nR10 and nOH) were profoundly chosen according to the structural features of the compounds and AChE active site. Furthermore, unlike some of the molecular descriptors usually used in QSAR models (Gurung et al., 2017, Wong et al., 2014), our variables are simple and easy to explain.

CRediT authorship contribution statement

Ante Miličević: Conceptualization, Methodology, Data curation, Writing – original draft, Writing – review & editing. Goran Šinko: Conceptualization, Data curation, Writing – original draft, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.