A quantitative structure-activity relationship study for structurally diverse HIV-1 protease inhibitors: contribution of conformational flexibility to inhibitory activity

Abstract

In this study, we investigated by linear regression model the SAR data of the 15 HIV-1 protease inhibitors possessing structurally diverse scaffolds. First, a regression model was developed only using the enzyme-inhibitor interaction energy as a term of the model, but did not provide a good correlation with the inhibitory activity (R² = 0.580 and Q² = 0.500). Then, we focused on the conformational flexibility of the inhibitors which may represent the diversity of the inhibitors, and added two conformational parameters into the model, respectively: the number of rotatable bonds of ligands (ΔSrot) and the distortion energy of ligands (ΔElig). The regression model by adding ΔElig successfully improved the quality of the model (R² = 0.771 and Q² = 0.713) while the model with ΔSrot was unsuccessful. The prediction for a training inhibitor by the ΔElig model also showed good agreement with experimental activity. These results suggest that the conformational flexibility of HIV-1 protease inhibitors directly contributes to the enzyme inhibition.

Keywords::

• Conformational flexibility
• distortion energy
• HIV-1 protease inhibitor
• linear regression model

Introduction

HIV-1 protease is one of the promising therapeutic targets of AIDS, and several protease inhibitor are widely and successfully being used for HIV/AIDS treatments 1-3. New inhibitors, however, are seamlessly being developed to overcome a drug-resistant virus that becomes an issue of the therapy 4-11.

The technology of structure-based design and/or computational chemistry has played a significant role in understanding the mechanism of HIV-1 protease inhibition and development of the protease inhibitors, e.g., identification of cyclic-urea inhibitors [12]. Among the previously reported QSAR studies on HIV-1 protease inhibitors, good correlation models of the inhibitory activity were obtained using the enzyme-inhibitor interaction energy. For instance, the Merck group analyzed 33 HIV-1 protease inhibitors and obtained good correlation models of the activity with the interaction energy using the MM2X/OPTIMOL method (R² = 0.7835 and Q² = 0.7551 as the best correlation among the models) [13]. Also, Gago and coworkers used the COMBINE methodology to analyze 49 inhibitors, and generated the highly predictive correlation model (R² = 0.91 and Q² = 0.81) [14].

Although the significant contribution of interaction energy to the potency of the inhibition has been demonstrated by these studies, the inhibitors used for the analyses were cognate series of inhibitors in terms of chemical structure. Therefore, one may argue about the limitation of designing new compounds having structurally diverse scaffolds using this model. In fact, the known protease inhibitors so far have a variety of chemical scaffolds, e.g., from linear to cyclic scaffolds. Recently, Karplus and coworkers successfully established the correlation model by the estimation of the absolute binding free energy of enzyme-ligand interaction for such structurally diverse inhibitors (R² = 0.83 and Q² = 0.71)[15]. This method is based on the conformational sampling by molecular dynamics and the electrostatic calculation by the Poisson or Poisson-Boltzmann equation. This model should be useful for designing a new compound with a new scaffold, but would be difficult to apply to medicinal-chemistry programs in non-computational labs.

In this study, we tried to apply a simple regression model for a set of structurally diverse protease inhibitors with the interaction energy alone initially, and then investigated the combination of other parameters to improve the regression model.

Material and methods

All forcefield-based computations were done by the MacroModel molecular modeling package with the BatchMin molecular mechanics engine [16]. The MMFF94s parameter was used as the forcefield [17]. The MMFF and AMBER94 [18] atomic charges were loaded for inhibitor and enzyme atoms, respectively. SGI O2 R5000 workstation was used for all computations. The crystallographic structures of HIV-1 protease complexed with inhibitors were obtained from the Protein Data Bank (PDB). The PDB entry IDs of the complexes and the inhibitory activities of the inhibitors are listed in Table I. The chemical structure of the inhibitors is shown in Figure 1.

Table I.  HIV-1 protease inhibitors used in this study.

Compounds	PDB ID	Ki (nM)	Number of rotatable bonds	Reference
A76889	1HVL	0.112	19	[34]
A76928	1HVK	0.011	19	[34]
A77003	1HVI	0.012	19	[34]
A78791	1HVJ	0.004	19	[34]
AG1343	1OHR	2.0	9	[35]
AHA001	1AJX	12.2	10	[36]
AHA006	1AJV	19.1	10	[36]
GR126045	1HTF	4.5	12	[37]
GR137615	1HTG	0.11	16	[37]
KNI272	1HPX	0.0055	14	[38]
L735524	1HSG	0.38	11	[39]
L738317	2BPV	21.2	13	[40]
SB203238	1HBV	430	14	[41]
SB206343	1HPS	0.6	16	[42]
U89360E	1GNO	20	19	[43]
JE2147	1KZK	0.041	9	[44]

Figure 1 Chemical structures of HIV-1 protease inhibitors used a) as a training set, and b) as a test.

Computational results and discussion

Preparation of enzyme-inhibitor complex structures

HIV-1 protease inhibitors used in this study were selected from SAR data published in various literature reports(Figure 1 and Table I). Although they all are classified into one class, peptidic-inhibitors, the SAR data should be appropriate for this study since the set of SAR data has diversity in terms of chemical structure, and the inhibitory activity (Ki) is widely dispersed from a single picomolar to subnanomolar level. Fifteen inhibitors were used for the QSAR analysis as a training set and one inhibitor (JE2147) was used for the prediction of its activity as a test.

To build a protease structure complexed with each inhibitor for computation, all water molecules were removed from the x-ray structure except for the water molecule Wat301, which is known to be commonly conserved in HIV-1 protease-inhibitor complexes, and to be hydrogen-bonded to the flip loops of the enzyme in a dimeric form [19]. For the complexes with cyclic-urea inhibitors, there is no water molecule equivalent to Wat301 since the carboxyl group of the inhibitors positionally replaces Wat301. Hydrogen atoms were then added to all heavy atoms of the enzyme and inhibitor, except for the aspartate residue at the 25th position of the enzyme (discussed below).

Determination of Asp protonation state

For the catalytically-significant aspartate residue at the 25th position, the side-chain atom of either Asp25 or Asp25′ in the dimer needs to be treated as a protonated form, i.e., one carboxyl and one carboxy anion groups in the side-chain atoms, according to known evidence [20]. To determine the protonation state of aspartate, two protonated states in a complex with each inhibitor were separately prepared, i.e., one complex has a protonation on Asp25 and another has a protonation on Asp25′. Both complexes were then energy-minimized by the conjugated gradient method with the 0.01 kJ/Å-mol gradient convergent criterion. The protonation state for a complex exhibiting lower conformational energy between the two minimized complexes was chosen as the protonation state of the aspartate for the complex (Table II). This assignment method has previously been used by Karplus's group [15]. Our results showed a good agreement with their assignment except for a few complexes including AG1343, L735524, and L738317.

Table II.  Energy (kcal/mol) of the complexes with protonation at Asp25 or Asp25′^a.

Compounds	Asp25	Asp25′	Compounds	Asp25	Asp25′
A76889	2082.01	2078.12	GR137615	2078.36	2083.89
A76928	2068.05	2072.16	KNI272	2085.45	2089.39
A77003	2090.63	2095.08	L735524	2080.32	2085.50
A78791	2073.64	2078.00	L738317	2086.01	2092.35
AG1343	2035.91	2033.57	SB203238	2100.22	2100.10
AHA001	2090.32	2088.66	SB206343	2100.48	2101.71
AHA006	2124.63	2122.15	U89360E	2027.08	2027.36
GR126045	2075.22	2074.35	JE2147	2059.17	2053.11

^a States with figures in italics denote the Asp protonation state to be chosen for calculation.

Calculation of enzyme-inhibitor interaction energy

For calculation of the enzyme-inhibitor interaction energy (ΔEint), the hydrogen-atom assigned complex structure was subjected to a short energy-minimization with the 0.1 kJ/Å-mol gradient convergent criterion in order to remove steric clash in the complex. For the minimization, the electrostatic energy was calculated by distance-dependant dielectric electrostatics with dielectric constant ϵ = 4r. The electrostatic non-bonded interaction was computed with a 12 Å cut-off distance and the van der Waals non-bonded interaction was computed with a 7 Å cut-off distance. Based on the short-minimized complex, ΔEint was calculated by the Equation (1). where E_complex is the conformational energy of the short-minimized complex; E_enzyme is the energy of the enzyme structure without the inhibitor; and E_inhibitor is the conformational energy of the inhibitor alone. The ΔEint values of each complex are listed in Table III.

Table III.  Calculated parameters for linear regression modeling.

Compounds	ΔEint (kcal/mol)	ΔSrot (kcal/mol)	ΔElig (kcal/mol)
A76889	− 60.91	5.915	5.982
A76928	− 65.47	5.915	7.223
A77003	− 65.42	5.915	5.417
A78791	− 61.70	5.915	4.441
AG1343	− 54.88	2.802	4.302
AHA001	− 49.75	3.113	1.766
AHA006	− 50.71	3.113	4.038
GR126045	− 46.21	3.736	2.297
GR137615	− 60.48	4.981	3.300
KNI272	− 54.25	4.358	4.026
L735524	− 60.90	3.424	5.157
L738317	− 54.62	4.047	7.321
SB203238	− 54.62	4.358	12.89
SB206343	− 56.19	4.981	3.490
U89360E	− 51.87	5.915	4.374
JE2147	− 59.45	2.802	5.226

Linear regression modeling with the interaction energy

To correlate the experimentally determined inhibitory activity (ΔGexp), linear regression equations were developed by the ΔEint value (Table IV). The squared correlation coefficient (R²) in the obtained regression model was not sufficiently high but the cross-validated correlation coefficient (Q²) was moderate (R² = 0.580 and Q² = 0.500). The regression model is shown as Model-1 in Table IV and calculated inhibitory activity (ΔGcalc) is listed in Table V. The quality of the Model-1 equation was not as good as the previously reported models, indicating that the SAR data in this study are not well explained by the interaction energy alone. A major difference in the SAR data between the previous and our studies is the diversity of the chemical structure of the inhibitors; the previous studies used cognate series of inhibitors, whereas the SAR data of this study were composed of structurally diverse inhibitors. Therefore, this suggests that a different regression model is needed for structurally diverse inhibitors.

Table IV.  Statistical results of linear regression models.

	ΔGexp	R^2	Q^2	SPRESS^a
Model-1	= 0.3139 × ΔEint +4.8756	0.580	0.500	1.594
Model-2	= 0.4407 × ΔEint + ΔSrot +7.4826	0.643	0.567	2.038
Model-3	= 0.3679 × ΔEint +0.4054 × ΔElig +5.8733	0.771	0.713	1.224

^a S_PRESS (kcal/mol) = the predictive residual error sum of squares.

Table V.  Experimental, calculated, and predicted free-energy differences.

Compounds	ΔGexp (kcal/mol)	ΔGcalc^a (kcal/mol)	Residual (kcal/mol)	ΔGcalc^b (kcal/mol)	Residual (kcal/mol)	ΔGcalc^c (kcal/mol)	Residual (kcal/mol)
A76889	− 14.16	− 14.25	− 0.086	− 13.45	0.712	− 14.11	0.048
A76928	− 15.59	− 15.68	− 0.082	− 15.46	0.137	− 15.29	0.308
A77003	− 15.54	− 15.66	− 0.120	− 15.43	0.106	− 16.00	− 0.459
A78791	− 16.22	− 14.49	1.726	− 13.79	2.424	− 15.03	1.192
AG1343	− 12.38	− 12.35	0.028	− 13.90	− 1.522	− 12.57	− 0.194
AHA001	− 11.26	− 10.74	0.519	− 11.33	− 0.069	− 11.71	− 0.453
AHA006	− 10.98	− 11.04	− 0.060	− 11.75	− 0.770	− 11.15	− 0.164
GR126045	− 9.82	− 9.63	0.189	− 9.15	0.672	− 10.20	− 0.377
GR137615	− 14.17	− 14.11	0.061	− 14.19	− 0.020	− 15.04	− 0.869
KNI272	− 16.02	− 12.15	3.870	− 12.07	3.957	− 12.45	3.570
L735524	− 13.21	− 14.24	− 1.030	− 15.93	− 2.721	− 14.44	− 1.231
L738317	− 10.53	− 12.27	− 1.741	− 12.54	− 2.013	− 11.25	− 0.725
SB203238	− 9.06	− 12.27	− 3.209	− 12.23	− 3.169	− 9.00	0.064
SB206343	− 13.12	− 12.76	0.359	− 12.30	0.821	− 13.39	− 0.263
U89360E	− 10.96	− 11.41	− 0.452	− 9.46	1.492	− 11.44	− 0.482
JE2147	− 14.78	− 13.79	0.996	− 15.91	− 1.134	− 13.88	0.902

^a Calculated using the Model-1 equation.

^b Calculated using the Model-2 equation.

^c Calculated using the Model-3 equation.

Linear regression modeling with conformational parameters

Structurally diverse compounds likely offer a variety of conformational flexibility, e.g., linear or fixed structure. In fact, the number of rotatable bonds in the SAR data of this study widely varied from 9 to 19 (Table I). Thus, to develop a good correlation model for the structurally diverse inhibitors, we focused on the conformational flexibility of the inhibitors and added conformational parameters into the Model-1 equation. So far, various conformational parameters to handle the conformational flexibility of ligands have been developed in virtual screening or QSAR modeling methods 21-27. Among them, in this study, two commonly used conformational parameters were employed to develop linear regression equations; one is a conformational entropic parameter derived from the number of rotatable bonds of ligands (ΔSrot), and the other is a conformational enthalpic parameter derived from the distortion energy of ligands (ΔElig).

ΔSrot represents the loss of conformational entropy, and was calculated by multiplying the number of rotatable bonds of a ligand (Nrot) by penalty coefficient (Crot) as shown in the Equation (2). where Crot = 0.3113 (kcal/mol/rotatable-bond) was used in this study. This value is used in Bohm's scoring function [23] and estimation of docking energy in AutoDock method [28]. In this study, the rotatable bonds of the inhibitors in the terminal functionalities and cyclic moieties were not counted. Since the ΔSrot value theoretically is positive or equal to zero, the following equation was applied for regression modeling so that the ΔSrot value is adopted as the penalty term in the free energy change equation.

The results of a linear regression equation with ΔSrot showed that the R² and Q² values were slightly improved to 0.643 and 0.567, respectively, from those of Model-1 (Model-2 in Table IV). This fairly small improvement could be due to the limitation of correct estimation of conformational flexibility only by the number of the rotatable bonds, e.g., difficulty to estimate flexibility of the cyclic moiety.

Next, the conformational enthalpic parameter, ΔElig, was used for the regression modeling. ΔElig was calculated by the subtraction of the conformational energy of a ligand in a local minimum (Elig_min) from the conformational energy of a ligand in the bound state (Elig_complex) as shown in Equation (3) (see Table III). The Elig_complex value used was the same as the E_inhibitor value defined in Equation (1). The Elig_min value was calculated by energy-minimizing the inhibitor structure, taken out from the complex, in free state.

As a result, a significant improvement of the correlation coefficients was achieved (R² = 0.771 and Q² = 0.713; Model-3 in Table IV). The quality of this regression model has become comparable to that of the previously reported models including the Karplus's model. The relevance of ΔElig in the binding process has been investigated by several computational studies 2129-32. The results obtained in this study also showed that the distortion energy of the inhibitors has a significant effect on the inhibitory activity for HIV-1 protease.

Prediction of inhibitory activity

The relevance of the obtained regression models, Model-1, Model-2, and Model-3, was examined by predicting the activity of an inhibitor, JE2147, which was excluded in the regression modeling (Figure 1 and Table I). The best result of the prediction was obtained from the Model-3 equation which has the highest Q² value and the lowest S_PRESS value among the regression models (Table V). The residual between the predicted and experimental activities for JE2147 was 0.902 kcal/mol.

Conclusion

In this study, the SAR data of HIV-1 protease inhibitors including structural diversity was investigated by linear regression models. We found that the linear regression model in a combination of the interaction energy with the distortion energy showed good correlation coefficients and predictive power while the interaction energy alone did not develop a good regression model. This suggests that the conformational flexibility of HIV-1 protease inhibitors directly contributes to the enzyme inhibition. In fact, in the field of medicinal chemistry, it is well known that the conformational flexibility of inhibitors significantly influences inhibitory activity, and the conformational fixation of molecules is one of the common tactics used to enhance biological activity [33].

In conclusion, the regression model deduced in this study has quality as good as the previous reported models, while it is composed only of two terms in the equation and heavy computation is not needed. Thus, this model can be useful for development of a new chemical scaffold for HIV-1 protease inhibitors.