Interpretable correlation descriptors for quantitative structure-activity relationships

Experimental

The data sets used in this study, 114 inhibitors of angiotensin converting enzyme [18] and 397 inhibitors of DHFR [17], have been widely used to investigate many QSAR methods [17]. Each data set contained experimentally determined activity (pIC₅₀) values for each molecule. TMACC descriptors were generated using the topological data of each molecule. All nonpolar hydrogen atoms were removed and their atomic value added to the heavy atom to which they were bonded. Polar hydrogen atoms were considered explicitly. Physicochemical properties were then assigned to each atom. Four parameters were used to represent these properties: Gasteiger partial charges [19], logS parameters [20], Crippen-Wildman molar refractivity parameters [21] and Crippen-Wildman partition coefficient (logP) parameters [21]. Gasteiger partial charges were calculated using the method of partial equalization of orbital electronegativity [19]. This procedure calculates atomic charges in σ-bonded and non-conjugated π-systems using only the topology of a molecule. LogS parameters were used to describe atomic contributions to aqueous solubility [20]. Crippen-Wildman molar refractivity (MR) was used as a measure of the steric effect, which is determined through classification of atoms based on neighbouring atoms [21]. Crippen-Wildman partition coefficients (logP) are assigned to each atom as a measure of atomic lipophilicity, determined in the same way as Crippen-Wildman molar refractivity.

The TMACC descriptor was derived from the autocorrelation descriptor [11]. The standard equation for calculating an autocorrelation descriptor, X _ac , (Equation 1), considers a property, p, and the topological distance, d, between atoms i and j:

(1)

TMACC descriptors are calculated as the product of the physicochemical properties as determined for every atom pair within a given molecule. Only the maximum value determined for any bond distance is used in the generation of the TMACC descriptors. All other values are neglected. As the TMACC descriptor incorporates both autocorrelation and cross-correlation, all possible combinations of physicochemical properties are considered. The equation for calculating a TMACC descriptor (X _TMACC ) (Equation 2) summarises this approach, involving two properties, p and q, for two atoms, i and j, separated by the topological distance, d:

(2)

All calculations were performed using Nottingham Cheminformatics Workbench (NCW), a package which provides the function of generating the TMACC descriptors and the TMACC interpretation. NCW is a Java-based application, which is suitable for all major operating systems. It builds upon the original TMACC code (available from our website http://comp.chem.nottingham.ac.uk/download/tmacc). The original software only generates TMACC descriptors; there is no facility to interpret your model. NCW allows the user to start with a set of molecules and perform a full analysis upon them. The popular machine-learning workbench Weka is included to provide PLS modeling, as well as an in-house implementation. The in-house PLS algorithm was written before Weka included one. The results of the PLS analysis are used to determine the atomic contribution of each atom. The interpretation is visualized by a colour scheme depicting activity contribution by atom. The user is able to view molecules individually or tabulated. NCW will be available as open source software for download from our website, http://comp.chem.nottingham.ac.uk/download/ncw. NCW utilizes Marvin for drawing and manipulating chemical structures: Marvin 5.2.2, 2009 http://www.chemaxon.com. All molecular graphics were created using YASARA ("Yet Another Scientific Artificial Reality Application," http://www.yasara.com).

Angiotensin converting enzyme (ACE) inhibition

LOO cross-validation of the PLS models generated for the ACE and DHFR data sets gave q² values of 0.70 and 0.53, respectively, consistent with those previously reported [9]. TMACC descriptors were generated for the ACE data set of 114 inhibitor molecules, previously used to assess the ability of a 3D QSAR method, comparative molecular field analysis (CoMFA) [18]. The data set contained a diverse range of structures and activities selected from literature [18]. ACE is a zinc metallopeptidase, which functions as a dipeptidyl carboxypeptidase, hydrolysing a range of oligopeptide substrates [22]. It acts to induce hypertension and is, consequently, a widely investigated target for antihypertensive drugs [23]. Several methods used to determine the ACE inhibitory activities (IC₅₀ values) involved the use of the substrate hippuryl-histidyl-leucine (HHL) [24, 25]. An early method for determining ACE inhibitor activity was to measure the rate of hippuric acid production from HHL catalyzed by ACE [26]. However, it was subsequently discovered that HHL is a C-terminal domain specific substrate of ACE [27, 28]. Consequently, the structure-activity relationship shown by this data set is likely to reflect that of C-terminal domain specific ACE inhibition, rather than general ACE inhibition.

For each molecule, a TMACC interpretation was generated, as described in the Methods section, leading to the labelling of each atom according to an activity banding (i.e., its contribution to activity). Some examples are shown in Figure 1. Based on the literature, several potentially important features (in the form of functional groups) were identified (Figure 2). Using the TMACC interpretation the extent of activity associated with each feature was accumulated for the entire ACE data set to allow determination of the structure-activity relationship modelled (Table 3).

	Activity
ACE inhibitor feature	Negative	Neutral	Positive
C-terminal carboxylate carbonyl	5	0	105
C-terminal carboxylate hydroxyl	0	5	105
Central carbonyl	13	3	96
Zinc binding carboxylate - carbonyl	3	3	22
Zinc binding carboxylate - hydroxyl	0	4	24
Zinc binding sulfhydryl sulfur	0	0	33
Zinc binding phosphinate phosphorus	0	0	22
Zinc binding phosphinate carbonyl	20	0	2
Zinc binding phosphinate hydroxyl	20	0	2
P1' methyl	4	2	27
P1' lysyl nitrogen	0	0	20

An essential feature of any ACE inhibitor is a zinc coordinating group. The catalytic zinc ion is coordinated by three highly conserved residues present in both somatic ACE (sACE) domains [29]. The important functional role of the zinc ion in the active site domains of ACE has led to the development of peptide based inhibitors, such as enalaprilat, with additional zinc-binding functional groups, including thiol, carboxylate and phosphinate groups. The importance of zinc binding functional groups in ACE inhibition has been demonstrated in crystal structure and structure-activity studies [22, 30, 31].

Zinc binding groups were frequently recognized by the TMACC interpretation as positive for activity. All sulfhydryl sulfur atoms located in the optimal position for zinc-binding were identified as positive for activity. Analysis of phosphinate zinc binding groups showed all phosphorus atoms were identified as positive for activity. However, phosphinyl oxygen atoms were mostly shown as negative for activity. In contrast, the interpretation most frequently identified both carboxylate zinc binding group oxygen atoms to be positive for activity. Although the results do not fully capture the correlation between the type of zinc-ligand and inhibitor activity observed in structure-activity studies, (phosphinate > carboxylate > sulfhydryl) [18], perhaps the negative activity attributed to the phosphinate oxygen atoms reflects its weak zinc-binding ability in comparison to the other zinc binding groups.

The central carbonyl group is a feature found in most ACE inhibitors. It forms two hydrogen bonds within both domains of ACE [22, 29]. Docking studies suggest this interaction is frequently present in ACE-inhibitor binding [32] and it has been identified in many ACE-inhibitor crystal structure complexes (Figure 3) [33]. Mutation of ⁵¹³His to alanine causes a 120,000-fold decrease in the binding of lisinopril to the C-domain of sACE [34]. This suggests the interaction of the conserved histidine residues with the carbonyl group of an inhibitor is important for ACE inhibition. The TMACC interpretation identified the central carbonyl as favourably contributing towards the activity (Table 3). The high frequency of positive activity shown for this feature by the TMACC interpretation is consistent with the aforementioned literature.

The crystal structures of testicular ACE (tACE) in complex with various inhibitors (Figure 3) show the intermolecular interactions responsible for ACE inhibition in tACE and correspondingly the C-terminal domain of sACE [22, 30, 33]. In contrast to most zinc protease inhibitors, which primarily rely on the strength of their zinc binding groups for activity, domain-specific ACE inhibitors utilize weak zinc binding groups and exploit both primed and unprimed sides of the active site in order to mimic peptide substrates, thereby achieving domain selective inhibition [35]. Domain-specific inhibition of ACE is important, as each domain possesses individual functions [36]. This discovery has developed the number of applications of ACE inhibitors, extending from treating hypertension to protecting stem cells during chemotherapy [37]. A recent study has also suggested ACE may be involved in many physiological processes other than blood pressure regulation [38].

A C-terminal carboxylate is found in many ACE inhibitors. This feature interacts with several conserved residues in both domains of sACE, hydrogen bonding with tyrosine and glutamine residues, and also forms an electrostatic interaction with a lysine residue (Figure 3) [32]. Both C-terminal carboxylate oxygen atoms were identified as positive by the TMACC interpretation (Table 3).

Despite the high level of conserved residues present in both domains of sACE, variation between the domains confers different substrate and inhibitor preferences. The presence of hydrophobic residues ³⁷⁹Val and ³⁸⁰Val in the S1' sub-site of the C-domain of sACE provides hydrophobic interactions between the sub-site and the P1' residue of inhibitor molecules, such as the P1' methyl group of captopril and enalaprilat [29]. The corresponding residues found in the N-terminal domain, ³⁵⁷Ser and ³⁵⁸Thr, provide a polar environment and, therefore, do not form similar hydrophobic interactions with the P1' residue of inhibitors [29, 33]. In the C-terminal domain the lysyl chain of lisinopril extends into the S1' sub-site and forms an electrostatic interaction with ¹⁶²Glu and a water-mediated interaction with ³⁷⁷Asp [22]. However, in the N-terminal domain the S1' sub-site makes fewer contacts with the lysyl chain of lisinopril (Figure 4). For example, ¹⁶²Glu (C-domain) is replaced by ¹⁴⁰Asp (N-domain), and due to the larger distance between the lysyl chain and this residue, no electrostatic interaction is observed at this location in the N-domain [29]. Additionally, ³⁷⁷Asp (C-domain) is replaced by ³⁵⁵Gln (N-domain), thereby abolishing the water-mediated interaction shown between the lysyl residue of lisinopril and the C-domain [29]. This evidence suggests that methyl and lysyl groups located in the P1' position of ACE inhibitors can form favourable interactions with the S1' sub-site of the C-terminal domain of sACE.

Interestingly, the TMACC interpretation identified all inhibitor P1' lysyl nitrogen atoms as favourably contributing towards activity (Table 3). The interpretation also identified inhibitor P1' methyl groups as positive for activity. Thus, the TMACC interpretation identified P1' groups important for C-domain specific ACE inhibition, as illustrated in Figure 1. This C-domain specific bias in the data set, reflected by the TMACC interpretation, has not been shown in previous QSAR investigations using this data set [17, 18].

Dihydrofoloate reductase (DHFR) inhibition

Due to the vital cellular function and ubiquitous expression of DHFR, inhibitors of this enzyme have been used clinically in the treatment of a range of diseases [39]. A variety of antifolates, which inhibit specific DHFR enzymes, have clinical application against cancer, malaria and many infectious diseases caused by bacteria, fungi and protozoa [40]. Natural folates contain a pteridine ring system, a p-aminobenzoic acid and a glutamate moiety (Figure 5). Classical antifolates are analogous to natural folates, possessing a glutamate residue, which allows polyglutamylation catalyzed by folylpolyglutamate synthetase [41]. Methotrexate (MTX) is a classical antifolate, which has been used clinically as an anticancer drug for over 50 years [42]. MTX is a potent inhibitor of DHFR from many species. However, it is restricted to anti-tumour applications, as a consequence of the dependence of classical antifolates on folate carrier-mediated transport, which is found only in mammalian cells [43]. Trimetrexate (TMQ) is a potent non-classical antifolate DHFR inhibitor, which is used in the treatment of Pneumocystis infections common in AIDS patients [44].

TMACC descriptors were generated for the DHFR data set of 397 molecules, previously studied [17]. The inhibitor activity values for the DHFR data set represent inhibitor potency against rat liver DHFR. Rat liver DHFR shows a high level of conservation of both primary sequence and active site residues and, therefore, is used as a mammalian standard for DHFR inhibition [45]. Following the strategy adopted for the ACE inhibitors, TMACC interpretations were developed for the DHFR inhibitors (Figure 6). Based on the literature, several DHFR inhibitor features were identified (Figure 7) and the association of each feature with activity was extracted from the TMACC interpretation and accumulated for the entire data set (Table 5).

	Activity
DHFR inhibitor feature	Negative	Neutral	Positive
2-amino nitrogen	155	81	161
4-amino nitrogen	110	61	212
N1-nitrogen	307	58	32
C8 Quinazoline	0	5	45
N8 Pyrido [2,3-d]pyrimidine	68	16	11
5-methyl	1	0	64
9-methyl	0	1	16
10-methyl	10	3	57
Benzylmethoxy - methyl	5	6	392

Table shows frequency of activity of DHFR inhibitor features as determined by the TMACC interpretation

An important feature for DHFR inhibition is a 2,4-diamino-substituted pyrimidine ring [49]. The 2-amino and 4-amino groups form hydrogen bonds with highly conserved residues, which orientate the pteridine ring of the inhibitor accordingly. The orientation of inhibitor binding differs from natural folate binding. Although both involve ³⁰Glu, the orientation of the inhibitor allows extensive hydrogen bonding with other DHFR residues, which is not possible in the folate binding orientation. Binding of inhibitors containing 2,4-diaminopyrimidine shows the 2-amino group forms a hydrogen bond with ³⁰Glu [50]. Additionally, the 4-amino group of the inhibitor forms strong hydrogen bonds with the carbonyl groups of ⁷Ile and ¹¹⁵Val [46, 50]. The protonated form of N1 forms an electrostatic interaction with ³⁰Glu [47, 51].

The N8 nitrogen of the pteridine ring forms a conserved network with a water molecule, which also hydrogen bonds with ³⁰Glu and ²⁴Trp [46, 52]. However, comparison of pyrido [2,3-d]pyrimidines and quinazolines suggests that the presence of the N8 hydrogen bond in inhibition by pyrido [2,3-d]pyrimidines may restrict the position of any bridge substituent, such as 9 or 10 methyl groups, preventing hydrophobic interaction of these features with the DHFR active site [53, 54]. Therefore, in general quinazoline DHFR inhibitors are more potent than pyrido [2,3-d] pyrimidine analogues.

The TMACC interpretation showed variable activity for the 2-amino group, with the activity approximately evenly distributed between negative and positive activity. However, interpretation of the 4-amino group showed a greater frequency of positive activity. The interpretation showed the N1 nitrogen as generally negative for activity. This classification could possibly result from the fact that the N1 was unprotonated within most molecules of the data set. However, studies suggest the protonated form of N1 forms a salt bridge with the DHFR active site [51, 55].

Comparison of the activity of C8 atoms in quinazoline inhibitors to N8 atoms in pyrido [2,3-d]pyrimidine inhibitors assigned by the TMACC interpretation showed the C8 atom of quinazoline inhibitors to be identified as positive for activity, whilst the N8 of pyrido [2,3-d]pyrimidines was most frequently negative for activity. This interpretation reflects experimental evidence, which indicates the presence of a C8 (quinazoline) is more potent for DHFR inhibition than an N8 (pyrido [2,3-d]pyrimidine) [56].

The presence of a 5-methyl group increases inhibitor potency by forming hydrophobic interactions with ¹¹⁵Val in human DHFR [53]. Methyl groups in the 9 and 10 positions of the bridge region also improve DHFR inhibitor potency, as the environment surrounding the bridge region is generally hydrophobic. Therefore, hydrophobic interactions with ²²Leu may be formed by these groups [47]. The benzyl ring featured in many antifolate inhibitors forms hydrophobic contacts with many residues within the DHFR active site [46, 51]. This feature is often substituted for other aromatic rings and methoxy groups to increase hydrophobic interactions.

The TMACC interpretation identified 5-methyl groups as positive for activity, consistent with experimental data, which suggests this feature is important for the inhibition of human DHFR. The interpretation also identified methyl groups in the 9 and 10 positions as positive for activity. This is consistent with the known hydrophobic interactions formed by these features within the DHFR active site [47]. The TMACC interpretation identified the methyl of methoxybenzyl groups substituted around the benzoic acid moiety of inhibitor molecules as positive for activity. This classification is supported by the hydrophobic interactions in which methoxy groups have been shown to participate within the DHFR active site [46, 51]. The TMACC interpretation of the DHFR data set identified many key structural features for DHFR inhibition. The analysis suggests that the hydrophobic groups investigated were more frequently identified as positive for activity than the hydrogen bonding groups investigated (Figure 6).