Docking, CoMFA and CoMSIA studies of a series of sulfonamides derivatives as carbonic anhydrase I inhibitors

Abstract

3D-QSAR methods, CoMFA region focusing (CoMFA-RF) and CoMSIA along with docking studies carried out for investigating 32 carbonic anhydrase I inhibitors. These inhibitors have been studied for the development of antiglaucoma, antitumor, antiobesity or anticonvulsant drugs. Docking analysis by GOLD provide conformations which have been realigned in CoMFA and CoMSIA models. Training set for the CoMFA–RF and CoMSIA models using 24 docked conformations gives q²_Loo values of 0.615 and 0.637 and r²_nv values of 0.701 and 0.713, respectively. The results of CoMFA-RF and CoMSIA with and without docked conformations were compared. The ability of prediction and robustness of the models were evaluated by test set, cross validation (leave-one-out and leave-ten-out), bootstrapping, and progressive scrambling approaches. The all-orientation search (AOS) was used to achieve the best orientation to minimize the effect of initial orientation of the structures. The docking results confirmed CoMFA and CoMSIA contour maps.

Keywords::

• Sulfonamides derivatives
• CoMFA
• CoMSIA
• docking

Introduction

Carbonic anhydrase(CA) (EC 4.2.1.1) is the first metalloenzyme identified in 1941 by Keilin and Man¹. These enzymes are belong to a family of zinc metalloenzymes, which catalyze the reversible hydration of carbon dioxide to bicarbonate ion (CO₂ + H₂O ⇔ HCO₃⁻ + H⁺)². So this is a basic reaction for respiration and transportation of CO₂ between metabolizing tissues and excretion sites³. CA enzymes have been classified in five distinct gene families: α-,β-,γ-, δ-, and ξ- that are present in prokaryotes and eukaryote^4,5. There are several cytosolic forms (CA I–III, VII), four membrane-bound isozymes (CA IV, IX, XII, and XIV), one mitochondrial form (CA V), and a secreted CA isozyme, CA VI⁶. The process of enzyme’s acting may be schematically represented by equations (1.1) and (1.2).

(1.1) (1.2)

The rate determining step is equation (1.2), that is, the proton transfer that reproduce the zinc hydroxide species of the enzyme⁷. The X-ray structures of CAs have shown that the zinc ion (Zn⁺²) within the active site is coordinated to three histidine residues and one water molecule (or hydroxide anion) and that the water molecule is replaced by a CA inhibitor as the fourth ligand of Zn⁺².⁸ One of the most prominent classes of CA inhibitors is aromatic/heterocyclic sulfonamides which have been studied for the development of antiglaucoma, antitumor, antiobesity or anticonvulsant drugs⁹.

Quantitative structure-activity relationship (QSAR) studies generally perform a vital role in drug discovery and design as ligand-based approaches¹⁰. Such approaches provide not only the reliable prediction of specific properties of new compounds, but also they help to clarify the possible molecular mechanism of the receptor-ligand interactions, in case that the experimental NMR or crystal structure of the target protein is not available¹⁰.

In recent years, a number of studies have been carried out for predicting new carbonic anhydrase inhibitors such as Clare and Supuran describe a QSAR based almost entirely on quantum theoretically calculated descriptors for a large and heterogeneous group of aromatic and heteroaromatic carbonic anhydrase inhibitors, using orbital energies, nodal angles, atomic charges, and some other intuitively appealing descriptors³. Gavernet and Coworkers synthesized sulfamides and sulfamates inhibitors and tested in their inhibition to carbonic anhydrase CAII activity, inhibition pattern interpreted by means of molecular modeling techniques⁴. Abdel-Hamid et al synthesized a new series of 1,3,4-thiadiazole-2-thione derivatives for the inhibition of three physiologically relevant carbonic anhydrase isozymes, the cytosolic human isozymes I and II, and the transmembrane, tumor-associated hCA IX and use docking studies of the tested compounds in order to predict the affinity and orientation of these compounds¹¹. Agrawal et al describe a QSAR study on a series of carbonic anhydrase inhibitors using a series of distance-based topological indices¹². Khadikar and coworkers describe a QSAR study on aromatic/heterocyclic sulfonamides containing 8-quinoline-sulfonyl carbonic anhydrase (CA) inhibitors which has been carried out topologically using first-order valence connectivity index (1vv)¹³. Melagraki and coworkers describe carbonic anhydrase II inhibitory with the method of linear quantitative structure–activity relationship for a series of para-substituted aromatic sulfonamides by using topological index methodologies¹⁴. Potter describes the docking of selected steroidal and non-steroidal estrone sulphatase inhibitors¹⁵. Quantum chemical QSAR expressions have been developed by Clare and Supuran for a heterogeneous group of 36 sulfonamides by using the ACE statistical technique¹⁶. Another Quantum chemical QSAR study which was performed by Clare and Supuran, have been developed on a group of 1,3,4-thiadiazole- and 1,3,4-thiadiazoline disulfonamides carbonic anhydrase inhibitors base on AM1 calculations¹⁷.

In this study we have investigated inhibition of a series of sulfonamides derivatives¹⁸ against carbonic anhydrase I by using 3D-QSAR models, comparative molecular field analysis (CoMFA)¹⁹ and comparative molecular similarity indices analysis (CoMSIA). In CoMFA, steric and electrostatic fields surrounding a set of aligned molecules in a grid box sampled and correlated with observed activities. In a similar method, CoMSIA, a probe atom is used to calculate similarity indices, at regularly spaced grid points for the aligned molecules. The molecular fields are calculated with potentials force field e.g. Lennard-Jones and Coulomb. While CoMSIA uses Gaussian-based similarity functions for molecular field calculations^19–21. Docking tools provide key structural features of binding of an inhibitor into the receptor and predicting bioactive conformers.

Experimental and methods

Data set

A series of 32 aromatic sulfonamide derivatives¹⁸ were investigated in this study Table 1. In vitro CAI inhibitory activities were changed into the corresponding pK_i (−log K_i) values. The set of inhibitors was divided into training and test sets. The test set compounds were selected by considering both the distribution of biological data and structural diversity of the molecules. At first the total set of CAI inhibitors (32 compounds) was divided into the training set (26 compounds) and test set (6 compounds), and 2 compounds (compounds no. 1 and 32) showed large residual values, then are identified as outliers and removed from the training set.

Table 1.  Structures and activities of aromatic sulfonamide derivatives.

Inhibitor	R	Ki(nM),hCA I
1	H	7.5
2	2-Me	107
3	3-Me	730
4	4-Me	104
5	2-F	621
6	3-F	116
7	4-F	108
8	2-Cl	640
9	3-Cl	311
10	4-Cl	112
11	2-Br	110
12*	3-Br	510
13*	4-Br	659
14*	3-OMe	342
15	F5	110
16	H	9
17	2-F	8.5
18	3-F	11.3
19	4-F	7.6
20	2-Cl	25.1
21	3-Cl	113
22	4-Cl	3.2
23	2-Br	43.4
24	3-Br	30.8
25	4-Br	12.3
26	2-Me	10.5
27*	3-Me	110
28*	4-Me	5.1
29	3-OMe	8.6
30	F5	9.7
31*		250
32		25

*Prediction set.

Molecular docking

The crystal structure of human carbonic anhydrase I (CA I) is 1AZM was taken from RCSB protein databank (http://www.pdb.org). This protein in PDB is not complexed with anyone of the understudy ligands, so in docking step its original ligand was removed and then ligands in our data set were docked in the active site of CA I one by one. Ligands preparation step was carried out in SYBYL 7.3 molecular modeling package (Tripos Inc., St. Louis, USA) running on a Red Hat Linux workstation 4.7. The resulting structures were imported into Discovery Studio 2.5 (Accelrys Inc, San Diego, CA, USA), and typed with CHARMm force field then partial charges were calculated by Momany-Rone option²². The resulting structures were minimized with Smart Minimizer which performs 1000 steps of steepest descent with a RMS gradient tolerance of 3, followed by conjugate gradient minimization²³. For preparation step of enzymes, all complexes were typed with CHARMm force field, hydrogen atoms were added, all water molecules were removed and pH of protein was adjusted to almost neutral, 7.4, using protein preparation protocol. All inhibitors were again minimized in-situ with Smart Minimizer option that is custom for in-situ ligand minimization and consists of some pre-defined minimization steps that have been pre-determined to work well for receptor ligand data²³. A 6.008 Å radius sphere was defined around the bounded ligands to confirm atoms of each ligand and the side-chains of the residues of the receptor within 6.008 Å from the centre of the binding site are free to move. Then bounded inhibitors were removed from the binding site. Other parameters were set by default protocol settings. GOLD program was used to dock inhibitors into receptor Supplementary Figure S1²⁴. Predicting how a small molecule will bind to a protein is complicated, and no method can guarantee success. The best approach is to measure as accurately as possible the reliability of the method (i.e., the chance that it will make a successful prediction in a given instance). For this reason, GOLD has been tested on a large number of complexes extracted from the Protein Data Bank. The overall conclusion of these tests was that the top-ranked GOLD solution was correct in 70–80% of cases²⁵. Different values of the genetic algorithm parameters may be used to control the balance between the speed of GOLD and the reliability of its predictions. GOLD will only produce reliable results if correct atom typing for both protein and ligand is used properly. In this method a site sphere around the ligand is needed and the radius was set to 8.954 Å in this approach.

CoMFA and CoMSIA

Two strategies were used to generate 3D-QSAR models, in first approach the structures of inhibitors were sketched in SYBYL molecular modeling package and partial atomic charges were calculated using the Gasteiger-Hückel method. The compounds were aligned to template molecule (compound 22) on a backbone which is common among all structures for minimizing the sum-of-squares deviation between reference backbone in each inhibitor and the corresponding core in the template. In CoMFA a probe sp³ carbon atom with +1 charge was employed to calculate steric and electrostatic interactions between the probe and structures. Electrostatic interactions calculated by using a Coulomb potential and van der Waals interaction using a Lennard-Jones potential function. Various column filtering values were also tested. CoMFA standard scaling applies the equal weight to data from each lattice point in any given field. Region focusing is an iterative procedure which refines a model by improving the weight for those lattice points which are most related to the model. This enhances the resolution and predictive capability (q²; cross validated r²) of a followed PLS analysis. Technically, this corresponds to rotate the model components during a high-order space²⁶. PLS region focusing is rationally equivalent to the GOLPE strategy and q²-GRS^27,28. In CoMSIA, standard settings (probe with charge +1, radius 1 Å and hydrophobicity +1, hydrogen-bond donating +1, hydrogen-bond accepting +1, grid spacing 2 Å) were used to calculate five different fields: steric, electrostatic, hydrophobic, hydrogen bond acceptor and donor. In 3D-QSAR, PLS analysis was used for modeling in which the independent variables were the CoMFA and CoMSIA fields, and −log K_i data values used as dependent variables.

In the second strategy bioactive conformers were docked and aligned inside the active site of protein by GOLD software, then all conformers from docking were realigned to the template molecule (compound 22) on a common backbone (Supplementary Figure S2) and other steps were performed like first strategy.

Because different space orientations of the molecular collective in the grid box have major effect on 3D-QSAR models, all-orientation search (AOS) was also carried out on initial orientations of aligned structures by the rotation procedure written in SYBYL programming language (SPL)²⁹.

Partial least squares (PLS) analysis

The CoMFA and CoMSIA descriptors were used as independent variables and pK_i values as dependent variables in partial least squares regression analysis. Leave one out(LOO) cross-validation method was used to obtain the optimal number of components (latent variables) in the subsequent analysis. The minimum-sigma (column filtering) was used to improve the signal-to-noise ratio in CoMFA and CoMSIA models. To get the highest correlation coefficient (r) and the lowest standard error in the LOO cross-validated predictions, the optimum number of principal components in the final non-cross-validated QSAR equations was determined. Analysis of CoMFA and CoMSIA results and the prediction of the models were performed by non-cross validation method according to SYBYL terminology.

Results and discussion

Docking results

Docking computations were employed to find the probable binding conformations of all carbonic anhydrase I inhibitors. To validate the docking reliability, root-mean-square distance (RMSD) value was calculated between bounded inhibitor and redocked ligand, which were 0.25 in this method. This value shows a high reliability of GOLD method to reproduce the known binding mode of these inhibitors. Docking results shows a hydrogen bond between SO₂NH₂ substituent and His119, Supplementary Figure S3. Another hydrogen bond was found between the proton of indole and HIS200, Supplementary Figure S3. According to docking results, there is an interaction between inhibitor molecules (in this case compound 22) and Zinc atom, by the oxygen of SO₂NH₂. Figure 1 shows three residues (His119, His96, and His94) and oxygen of SO₂NH₂ which have interaction with Zinc atom. MOLCAD surfaces were calculated for the most active compound (M22) to demonstrate electrostatic potential (Supplementary Figure S4a), lipophilic potential (Supplementary Figure S4b) and hydrogen bonding (Supplementary Figure S4c) interactions.

Figure 1.  Interaction among Zinc, protein residues and ligand.

CoMFA and CoMSIA results

PLS analysis of the compounds in training set showed CoMFA-region focusing (CoMFA-RF) QSAR model (grid spacing = 1) with a good q² value of 0.627 (1 components) that is better than common CoMFA. The optimal number of components evaluated by selecting the highest q² value corresponds to lowest S_press value. The non cross-validated PLS analysis results in a high conventional showed r² of 0.706, F = 52.803 a low standard error of estimation (SEE) 0.384, with a column filtering of 3.0.

The CoMSIA analysis was done at a grid spacing 2 Å, and the effect of column filtering was checked with the combination of five fields. The CoMSIA method defines hydrophobic and hydrogen bond donor and acceptor descriptors in addition to the steric and electrostatic fields in CoMFA. The combination of fields was systematically changed to select the optimal results.

By using the combination of five fields, q² of 0.636 was obtained with 1 component at a column filtering of 2.5 kcal/mol, The non cross-validated analysis results showed r² of 0.713, F = 54.597, and SEE = 0.380. The q² values of each independent field of steric, electrostatic, hydrophobic, hydrogen bond donor and hydrogen bond acceptor are 0.648, 0.608, 0.619, 0.593, 0.645 respectively. The contribution of all of these five fields has significant effect on constructed model. The values of experimental and predicted activities are depicted in Table 2.

Table 2.  The experimental K_i values, predicted K_i values (Pred.) and the residuals (Res.) of the training and test set compounds based on realignment method.

Compound no.	Experimental	CoMFA-RF		CoMSIA
		Pred.	Res.	Pred.	Res.
1	8.125	6.697	1.428	6.694	1.431
2	6.971	6.737	0.234	6.687	0.284
3	6.137	6.737	−0.600	6.702	−0.565
4	6.983	6.736	0.247	6.694	0.289
5	6.207	6.706	−0.499	6.673	−0.466
6	6.936	6.681	0.255	6.729	0.207
7	6.967	6.678	0.289	6.725	0.242
8	6.194	6.666	−0.472	6.719	−0.525
9	6.507	6.683	−0.176	6.726	−0.219
10	6.951	6.717	0.234	6.752	0.199
11	6.959	6.726	0.233	6.746	0.213
12*	6.292	6.734	−0.442	6.758	−0.466
13*	6.181	6.730	−0.549	6.755	−0.574
14*	6.466	6.753	−0.287	6.736	−0.270
15	6.959	6.74	0.219	6.731	0.228
16	8.046	7.901	0.145	7.903	0.143
17	8.071	7.911	0.160	7.937	0.134
18	7.947	7.909	0.038	7.909	0.038
19	8.119	7.906	0.213	7.912	0.207
20	7.600	7.665	−0.065	7.701	−0.101
21	6.947	7.667	−0.720	7.692	−0.745
22	8.495	7.685	0.810	7.695	0.800
23	7.363	7.936	−0.573	7.951	−0.588
24	7.511	7.925	−0.414	7.911	−0.400
25	7.910	7.934	−0.024	7.913	−0.003
26	7.979	7.928	0.051	7.943	0.036
27*	6.959	7.902	−0.943	7.876	−0.917
28*	8.292	7.916	0.376	7.901	0.391
29	8.066	7.751	0.315	7.586	0.480
30	8.013	7.914	0.099	7.904	0.109
31*	6.602	6.942	−0.340	6.947	−0.345
32	7.602	6.990	0.612	7.063	0.539

*Prediction set

Beside the 3D-QSAR studies described above that constructed and aligned conformers in SYBYL program, all conformers were docked in CAI protein by using GOLD method in discovery studio program and then realigned to the template molecule (compound 22 in this case) in SYBYL program by using rigid alignment method. A CoMFA-RF 3D-QSAR model with a q² value of 0.615 with 1 PLS components obtained. Also, CoMSIA 3D-QSAR study was performed on the compounds in the realignment method. CoMSIA PLS analysis determined q² value of 0.637 with 1 PLS components. Table 3 shows results of models constructed by these two alignment methods.

Table 3.  Summary of the statistical results for the constructed models.^a

Statistical parameters	Rigid alignment		Realignment
	CoMFA-RF (cf = 3, vars = 1728)	CoMSIA (cf = 2.5, vars = 1728)	CoMFA-RF (cf = 3, vars = 1560)	CoMSIA (cf = 3, vars = 1560)
q^2	0.627	0.636	0.615	0.637
SEP	0.433	0.427	0.440	0.427
r^2ncv	0.706	0.713	0.701	0.713
SEE	0.384	0.380	0.387	0.379
Fratio	52.803	54.597	51.664	54.766
r^2pred	0.744	0.706	0.706	0.702
Component	1	1	1	1

^aThe q² values for all PC’s from 1 to 6 are the same and the PRESS value for the first PC’s is minimum.

Validation of the 3D-QSAR models

The external set of 6 compounds was used to confirm predictive ability of the models. The r²_pred from CoMFA-RF and CoMSIA models without using docking conformer were found to be 0.744 and 0.706, and with using docking conformer were 0.706 and 0.702 respectively, which show models have acceptable predictability. To evaluate the statistical confidence limits of the derived models, bootstrapping¹⁹ analysis was carried out with 100 runs. Bootstrapping subsumes the generation of many new datasets from the original dataset after randomly choosing samples from that. A r²_bs (average correlation coefficient for bootstrapping) of 0.735 ± 0.071 and a SEE_bs (average standard error of estimate for bootstrapping) of 0.353 ± 0.182 for CoMFA-RF model with using docking conformer, and r²_bs of 0.738 ± 0.072, and a SEE_bs of 0.358 ± 0.144 for CoMSIA model generated by docking conformer, suggested a good internal consistency and the absence of systematic errors of the models. To evaluate the sensitivity of the optimized CoMFA-RF and CoMSIA models to chance correlations, the leave-one-out (LOO), leave 10-out cross-validation and progressive scrambling analyses were performed³⁰. The q² values of leave 10-out for CoMFA-RF and CoMSIA models without using docking conformer were 0.630 and 0.633 and with using docking conformer were 0.611 and 0.635 respectively. In the progressive scrambling approach, small random perturbations are introduced into a data set and the statistical results, the perturbation prediction (q²), the calculated cross-validated standard error of prediction (cSDEP) as the function of the correlation coefficient between the true values (y) of the dependent variables and the perturbed values (y’) of the dependent variables, and the slope of q² (cross validated correlation coefficient) with respect correlation of the original dependent variables against the perturbed dependent variables (d_q^2′/dr²_yy’), for CoMFA-RF and CoMSIA models are summarized in Table 4.

Table 4.  Model progressive scrambling for CoMFA-RF and CoMSIA models based on realignment method.

Model	q^2	cSDEP	dq^2′/dr^2yy’
CoMFA-RF	0.475	0.478	0.900
CoMSIA	0.490	0.525	0.973

CoMFA contour maps analysis

In the CoMFA steric contour maps, green contours show regions that are sterically favorable and yellow contours show sterically unfavorable regions, these contours represent 80% and 20% contributions, respectively. Similarly blue contours show the regions that are electropositive charge favorable and red contours show electronegative charge favorable regions that their contributions in the CoMFA electrostatic field represent 80% and 20% respectively.

In the first series (compound 1–15) there is a steric unfavorable yellow contour cover R₃ and R₂ substituents indicate that bulky groups in these regions will decrease inhibitory activity Supplementary Figure S5. It can be demonstrated by comparing the structures and activities of compound 3 (R₃ = Me, K_i = 730) with 6 (R₃ = F, K_i = 116) and compound 5 (R₂ = F, K_i = 621) with 8 (R₂ = Cl, K_i = 640). In the second series (compound 16–32) near R₂ and R₃ substituents there are steric favorable green contours indicate that in these positions bulky groups increase inhibitory activity Figure 2a. For example compound 29 (R₃ = OMe, K_i = 8.6) is more active than compound 27 (R₃ = Me, K_i = 110).

Figure 2.  CoMFA contour map displaying steric (a) and electrostatic (b) in combination with compound 27 and 22 respectively based on realignment method.

In CoMFA electrostatic contour maps for second series, there is a big red contour near the substituent of pyridine ring (Figure 2b) that shows electronegative groups in this situation increase the inhibitory activity, so that activity of second series that has N⁺ in pyridine ring is more than first series which has not N⁺. For example compound 22 (R₄ = Cl, K_i = 3.2) is more active than compound 10 (R₄ = Cl, K_i = 112).

CoMSIA contour map analysis

The steric field distribution of CoMSIA model is shown in Figure 3a. By comparison of Figure 3a with Figure 2a, it can be seen that CoMSIA steric contour maps is similar to CoMFA steric contour maps. In this field for second series there is a big green contour covers pyridine ring indicates that bulky groups in this region increase activity. It can be demonstrated by comparing the activity of second series with first series, for example compound 17 (R₂ = F, K_i = 8.5) which is in second series is more active than compound 5 (R₂ = F, K_i = 621) in first series.

Figure 3.  CoMSIA contour map displaying steric (a), electrostatic (b), hydrophobic (c), hydrogen bond donor (d), hydrogen bond acceptor (f) in combination with 17, 17, 16, 16 and 23 respectively based on realignment method.

In the CoMSIA electrostatic contour maps in addition to the contours that were shown in CoMFA, for second series there is a big red contour which covers R₂ substituent indicates that electronegative groups in this position increase inhibitory activity, Figure 3b. It can be demonstrated that compound 17(R₂ = F, K_i = 8.5) is more active than compounds 20(R₂ = Cl, K_i = 25.1) and 23(R₂ = Br, K_i = 45.5).

The CoMSIA hydrophobic contour maps are shown by Figure 3c. For this field the yellow (hydrophobic favorable) contours and white (hydrophobic unfavorable) contours represent 80% and 20% contributions, respectively. For second series compounds (16–32) there is one big yellow contour covers the phenyl ring which is substitutated in 3-position of indole indicates that hydrophobic groups in this area are preferred for CAI inhibitory activity. According to docking results in this area there are some interaction with hydrophobic residues of receptor, such as LEU131, PHE91 and GLY92 that confirm CoMSIA’s hydrophobic contour Supplementary Figure S6. This is a sensible reason why compound 16 (R = H, K_i = 9.00) has higher activity than compounds 21 and 27 (R = Cl, K_i = 113 and R = Me, K_i = 110).

Figure 3d depicts H-bond donor field distribution of CoMSIA model which is represented by cyan and purple contours. Cyan and purple contours indicate regions where hydrogen bond donor substituents on ligand are favored and disfavored respectively. There is one purple contour cover the−CO group in CONH₂ indicate that hydrogen bond donor groups in the receptor enhance the activity of inhibitor, Supplementary Figure S7. Based on hydrogen bond acceptor field in Figure 3f magenta contours show regions where hydrogen bond acceptor groups are favored and red contours indicate regions where hydrogen bond acceptor groups are unfavorable for increasing activity. There is one magenta contour near carbonyl substituent suggest hydrogen bond acceptor in this region is favored. Complementary of this magenta contour in the receptor, is hydrogen bond donor group (His 200). Also presence of electrostatic red contour at this position confirms the magenta contour.

Conclusions

In this study, molecular docking and 3D-QSAR studies were performed on a series of sulfonamide CAI inhibitors. 3D-QSAR models were generated by two different methods, by using docked ligands as a bioactive conformation and without docked ligands. The comparison of these models indicates that CoMFA and CoMSIA models which were generated by docked conformers have more reliable contour maps. The CoMFA and CoMSIA analyses have provided discerning key structural features which affect inhibitory activity of these inhibitors. We identified and confirmed the residues that play key role in the hydrogen bond donors and hydrogen bond acceptors in addition to the hydrophobic, steric and electrostatic interactions. The significant statistical parameters and appropriate predictive ability of the generated models indicate that these models can help to rational design of novel carbonic anhydrase inhibitors.

Declaration of interest

The authors report no conflicts of interest.