a
Department of Chemical Pathology, Faculty of Health Sciences, University of Pretoria
and National Health Laboratory Service Tshwane Academic Division, Pretoria, South
Africa.
b
Division of Chemical Pathology, University of Cape Town, , South Africa.
Correspondence should be addressed to T.S. Pillay, Department of Chemical Pathology,
Faculty of Health Sciences, University of Pretoria, Private Bag X323, Arcadia, Pretoria,
0007
Email: tspillay@gmail.com
Phone: +27-12-319-2114
Fax: +27-12-3283600
Structural requirements for potential HIV-integrase inhibitors
identified using pharmacophore-based virtual screening and
molecular dynamics studies†
Md Ataul Islam
a
, Tahir S. Pillay
*a,b
Acquired immunodeficiency syndrome (AIDS) is a life-threatening disease which is a collection of symptoms and infections
caused by a retrovirus, human immunodeficiency virus (HIV). There is currently no curative treatment and therapy is
reliant on the use of existing anti-retroviral drugs. Pharmacoinformatics approaches have already proven their pivotal role
in the pharmaceutical industry for lead identification and optimization. In the current study, we analysed the binding
preferences and inhibitory activity of HIV-integrase inhibitors using pharmacoinformatics. A set of 30 compounds were
selected as the training set of a total 540 molecules for pharmacophore model generation. The final model was validated
by statistical parameters and further used for virtual screening. The best mapped model (R = 0.940, rmsd = 2.847, Q
2
=
0.912, se = 0.498, R
2
pred
= 0.847 and r
2
m (test)
= 0.636) explained that two hydrogen bond acceptor and one aromatic ring
features were crucial for the inhibition of HIV-integrase. From virtual screening, initial hits were sorted using a number of
parameters and finally two compounds were proposed as promising HIV-integrase inhibitors. Drug-likeness properties of
the final screened compounds were compared to FDA approved HIV-integrase inhibitors. HIV-integrase structure in
complex with the most active and final screened compounds were subjected to 50ns molecular dynamics (MD) simulation
studies to check comparative stability of the complexes. The study suggested that the screened compounds might be
promising HIV-integrase inhibitors. The new chemical entities obtained from the NCI database will be subjected to
experimental studies to confirm potential inhibition of HIV integrase.
Introduction
Human immunodeficiency virus (HIV) is the aetiological agent of
acquired immunodeficiency syndrome (AIDS) which destroys the
immune system of the body leaving the victim vulnerable to
infections, malignancies and neurological disorder. Owing to its
rapid spread it has become a serious global threat and there is no
curative treatment for this fatal disease. According to statistics by
World Health Organization (WHO), a total of 37.20 million people
are living with AIDS and 1.70 million people died in 2013 alone.
Currently, there are 3 categories of therapeutic anti-HIV drugs
based on their inhibitory mechanisms
1
and these include
nucleoside reverse transcriptase inhibitors (NRTIs)
2
, non-nucleoside
reverse transcriptase inhibitors (NNRTIs)
3
, and protease inhibitors
(PIs)
4, 5
. To date the highly active antiretroviral therapy (HAART)
6
which is combined therapy using the above classes of inhibitors is
widely used for patients with advance infection but has failed to
eradicate the virus. HAART is intended to slow down viral
replication and lower the patient’s total burden of HIV infection,
but this treatment is not entirely cost effective and is often out of
reach of people worldwide. The genome of the HIV encodes for
three enzymes viz. the protease, reverse transcriptase and
integrase. The integrase is a 32 kDa enzyme made of three
functional domains included an N-terminal domain, catalytic core
domain and a less conserved C-terminal domain
7, 8
. The HIV
1
integrase has no equivalent counterpart or sequence homologue in
the human host cell and consequently it can be considered as an
attractive drug target
9
. After the approval of Raltegravir
10
as anti-
HIV drug several integrase inhibitors emerged as promising class of
therapeutics for the treatment of AIDS. Raltegravir and several
other inhibitors have been identified to possess anti-HIV integrase
activity but these have adverse effects on prolonged use and the
development of drug resistance drives the need to explore new
novel and potential chemical scaffolds for the treatment of AIDS.
The computational methods in drug discovery collectively termed
pharmacoinformatics, includes structure activity relationship (SAR),
pharmacophore, virtual screening and molecular docking and these
have proven their pivotal role in the pharmaceutical industry for
lead identification and optimization. Several research groups
worldwide identified integrase inhibitors
11-17
using
pharmacoinformatics approaches for potential application for HIV
therapy. Consistent with the objective of developing new potent
and less toxic integrase inhibitors the current research explores the
binding preferences of the inhibitory molecules of HIV integrase in
terms of space modelling study and virtual screening along with
molecular docking and molecular dynamics.
A pharmacophore model is an collection of steric and electronic
features and provides an intuitive way of depicting and
understanding the binding properties of small molecules along with
an explanation of optimum supra-molecular interactions with a
precise biological target, to activate (or block) its biological
response
18, 19
. It can also be defined that the pharmacophore idea is
based on the kinds of interaction observed in molecular
appreciation, i.e., hydrophobic, hydrogen bonding, and charge
interaction. For the HIV integrase inhibitors the hydrogen bond
acceptor (HBA) and donor (HBD), hydrophobic (H) and aromatic ring
(RA) pharmacophore features were found to be the important
functional features associated with the selectivity and potency. The
pharmacophore models are widely used in the field of drug
discovery by providing valuable information to study SAR and
reveals the mechanism of ligand-target relationship by deducing the
nature of functional groups and non-covalent bonding patterns
20
. It
can also be used in virtual screening to identify potential molecules,
predict the activity of the newly synthesized compound before
animal experiment; or understand the possible mechanism of
action
21, 22
. In the current study, an attempt was made to identify
the pharmacophore hypothesis using the HypoGen module
23
based
on key chemical features of HIV-integrase inhibitors with inhibition
constant covering a satisfactory wide range of magnitude. The
model was validated using several statistical approaches including
Fischer’s randomization and test set prediction. The validated
model was utilized for the virtual screening to select the virtual hits
from structural database. The molecular docking study was also
performed to elucidate the binding interactions and preferred
orientation of proposed potential molecules. The potential of the
work is displayed by the identification of two potential lead
molecules as integrase inhibitors. Finally, the molecular dynamics
study was performed to analyse stability and precise binding
interactions of the screened molecules inside the receptor cavity of
HIV integrase.
Materials and methods
At present, several popular commercial and freeware packages are
used for ligand-based method to derive 3D pharmacophore models
and also help in estimation of biological activities. Here we used
Discovery Studio 4.0 (DS)
24
for the 3D QSAR pharmacophore, virtual
screening and molecular docking studies. This is commercially
available software containing several module packages and widely
used in pharmacoinformatics drug discovery
25-28
. The 3D QSAR
Pharmacophore Generation module takes input of structure and
activity data for a set of potential HIV-integrase ligand to create
hypotheses. Two modules, HypoGen and HipHop are used for
ligand-based pharmacophore modelling. The HypoGen allows
identification of hypotheses that are common to the ‘active’
molecules of training set but absent in the ‘inactive’ molecules,
whilst HipHop identifies hypotheses present both in ‘active’ and
‘inactive’ compounds. In the present work the HypoGen module
was used to generate the hypotheses.
Dataset
1437 compounds belong to a collection of HIV-1 integrase inhibitors
were downloaded from BindingDB (http://www.bindingdb.org/)
with data on inhibitory activity (IC
50
). Duplicate and compounds
without definite activity values were deleted and finally 540
compounds considered as whole dataset for the study. Training and
test set compounds were separated from whole dataset for
pharmacophore model generation and validation of generated
model respectively. The molecules of the dataset have a wide range
2
of IC
50
, from 2.000 to 1000000.000 nM. For simplicity the whole
dataset was divided into three sets on the basis of inhibitory
activities values; highly active (IC
50
< 100.000 nM, +++), moderately
active (100.000 ≤ IC
50
< 1000.000 nM, ++) and least active (IC
50
≥
1000.000 nM, +). For selection of the training set for
pharmacophore model generation in DS basic guidelines laid down
by Li et al.
29
were followed. The guidelines reported as a) molecules
should be selected to provide clear and brief information with
structure features and range of activity, b) at least 16 diverse
molecules for training set should be considered to ensure the
statistical significance and avoid chance correlation, c) the training
set must include the most and the least active molecules and d) the
biological activity data of the molecules should have spanned at
least 4 orders of magnitude. Following the above guidelines four
training sets were (Set 1, Set 2, Set 3 and Set 4) generated
containing 30 compounds each. It was also kept in mind that no
compounds were common in any two training sets except for the
most active and least active molecules. The remaining 410
molecules were considered as test set molecules for each set and
used for assessing the performance of pharmacophore model. The
2D/3D visualizer
24
of DS was used to generate three-dimensional
coordinates of the compounds. For each compound, the
coordinates were corrected, atoms were typed and energy was
minimized using the modified CHARMm force field
30, 31
. The several
packages of DS were used for pharmacophore, virtual screening and
molecular docking studies.
Pharmacophore model generation
In order to generate the pharmacophore space model the 3D QSAR
Pharmacophore Model Generation module of DS was used.
Conformations of the training set molecules were generated by Cat-
Conf program of the DS software package. Out of BEST/FAST, the
BEST method was considered to obtain multiple acceptable
conformations which provides complete and enhanced coverage of
conformational space with help of rigorous energy minimization
and optimizing the conformations by the poling algorithm
32
. In the
BEST algorithm, the chemical features are arranged in space instead
of simply the arrangement of atoms
33
. For prediction of the
favourable features for the highly active compounds of the dataset
the Feature mapping was considered. Mapped features were given
as input features for pharmacophore model generation. Using the
conformer along with chemical features the modules operates in
two modes such as HipHop and HypoGen. The HipHop approach
generates the pharmacophore models by using active compounds
only, while the HypoGen approach considered both active and
inactive compounds in order to find out a hypothesis which is
common in the active molecules and absent in the inactive
compounds
33
. Top ten hypotheses are generated by the HypoGen
with consideration of the training set, conformational models and
chemical features through three steps: constructive, subtractive
and optimization
34
. In the first step, hypotheses are generated that
are common in the most active compounds; in subtractive phase,
inactive compounds are removed from those that fit the
hypotheses. In final step, the remaining hypotheses improve the
score with help of small perturbations
33, 35
. The best hypothesis was
selected based on the best correlation coefficient (R), low root
mean square deviation (rmsd), cost function analysis and good
predictive ability.
All four sets were used to develop the pharmacophore models and
statistical parameters were calculated based on training & test sets
molecules. Statistical results are depicted in Table S2 in
Supplementary file. It was observed that Set 1 gives the better
statistical results compared to Sets 2 – 4. Hence the Set 1 (n
tr
= 30,
Fig. 1) was considered as the training set in the current study. In the
remaining section “training set” is explained as Set 1 compounds
and “test set” as test compounds of Set 1. Training set molecules in
SMILES format with observed and estimated activity along with fit
values of Set 2, Set3 and Set 4 are depicted in Tables S2, S3 and S4
respectively in Supplementary data. The information concerning the
structure and the biological activity of test set compounds of Set 1
is provided in Table S5 in the supplementary information, while all
the data regarding the training set (Set 1) molecules are reported in
Fig. 1. The activity value distribution of the training set molecules of
Set 1 is given in Figure S1 (Supplementary file).
3
Fig. 1 2D chemical structures of the training set compounds and the
activity values (IC
50
) are given in the parentheses.
Validation of pharmacophore model
In order to check the predictivity and applicability as well as
robustness of the pharmacoinformatics model, the pharmacophore
hypotheses developed were validated by five different methods, (1)
internal validation, (2) cost function analysis, (3) Fischer’s
randomization test, (4) test set prediction and (5) decoy set.
Internal validation
Leave-one out (LOO) cross-validation is one of the important
internally validation protocol of the selected model, in which one
compound was randomly deleted from training set in each cycle
and model redeveloped using the rest of the compounds with the
same parameters used in original model. The activity of the deleted
compounds was calculated based on the newly developed model.
The above procedure applied for all molecules of the training set
and predicted activity recorded. Two important statistical
parameters, the LOO cross-validated correlation coefficient (Q
2
) and
error of estimation (se) were calculated based upon predicted
activity of training compounds. It is reported that high Q
2
(>0.5) and
low se explained better predictive ability
36
.
Another parameter, the modified r
2
(r
2
m(LOO)
) developed by Roy et
al.
37, 38
was calculated to confirm the good predictive ability of the
training set molecules. This parameter measure the degree of
deviation of the predicted activity from the observed ones. It was
stated that model may be considered with r
2
m(LOO)
>0.5.
Cost function analysis
To select the final pharmacophore model several statistical
parameters were employed at the time of hypothesis generation,
these included spacing, uncertainty, and weight variation. The
spacing represents the minimum inter-features distance which may
be permissible in the final hypothesis, while the weight variation
reflects the level of magnitude explored by the hypothesis in which
every feature implies some degree of magnitude of the biological
activity of the compound. The default values of spacing and weight
variation are 300 and 0.3 respectively but some cases it varies from
400 to 100 and 1 to 2 respectively. The uncertainty returns the
error of prediction which signifies the standard deviation of the
error cost. The default value of uncertainty parameter is 3 but it
may vary from 1.5 to 4.0 depending upon the nature of work. The
cost-function analysis is an important aspect for selection of final
model which minimized three terms, viz., weight cost, error cost,
and configuration cost. The weight cost is directly proportional to
the deviation of weight variation from its input value. The error cost
is the deviation between the predicted activity and experimentally
determined activity of the molecule in the training set. A fixed cost
is determined by the complexity of the hypothesis space. The
configuration cost implies entropy of hypothesis space and it is
reported that value should have <17 for a good pharmacophore
model. HypoGen module also calculates the null hypothesis which is
the assumption that there is no relationship in the data, and the
experimental activities are distributed about their mean. It is
N
N
OH
O
HN
F
O
N
O
H1 (60.000 nM)
N
N
OH
O
HN
F
O
N
H2 (62.000 nM)
N
OH
O
HN
O
F
N
N
S
O
O
N
H3 (5.000 nM)
N
N
OH
O
HN
F
O
N
H
N
S
O
O
H4 (30.000 nM)
N
N
OH
O
HN
F
O
N
N
O
N
H
H5 (52.000 nM)
N
N
OH
O
HN
O
N
O
F
Cl
H6 (28.000 nM)
N
OH
O
HN
O
F
N
N
H7 (4.000 nM)
N
N
OH
O
HN
F
O
N
N
S
O
O
H8 (11.000 nM)
N
OH
O
HN
O
F
N
N
O
H9 (5.000 nM)
N
OH
O
HN
O
F
N
N
O
N
N
H10 (3.000 nM)
N
N
OH
O
HN
O
N
O
F
Br
H11 (44.000 nM)
N
OH
O
HN
O
F
N
N
S
O
O
N
H12 (10.000 nM)
N
O
OH
O
H
N
F
O
N
N
H13 (2.000)
N
N
OH
N
N
O
F
S
HN
O
O
H14 (12.000 nM)
N
N
OH
N
N
N
O
F
N
O
O
N
H15 (3.000 nM)
N
N
OH
N
N
O
F
O
H
2
N
H16 (6.000 nM)
N
N
OH
O
N
H
F
H17 (12.000 nM)
N
N
OH
N
N
O
F
N
N
H18 (20.000 nM)
N
N
OH
N
N
N
O
F
S
O
O
H19 (8.000 nM)
N
N
OH
N
N
N
H
H20 (2400.000 nM)
O
O
O
OH
F
H21 (250.000 nM)
O
O
O
OH
O
H22 (140.000 nM)
O
O
O
O
O
O
O
O
HO
OH
H23 (1900.000 nM)
N
H
H
N
S
S
O
O
O
O
O
O
H24 (97000.000 nM)
N
N
HN
H
N
N
O
O
O
HO
F
H25 (12.800 nM)
HN
NH
OH
O
HO
O
O
O
H26 (28000.000 nM)
N
N
H
N
O
HO
F
N
O
H27 (25.000 nM)
NH
NH
HO
OH
H28 (10000.000 nM)
H29 (1000000.000 nM)
N
N
OH
O
O
O
F
Cl
O
O
H30 (15600.000 nM)
N
O
O
O O
Cl
F
4
illustrated that the higher (>60) cost difference (∆cost = null cost -
total cost) indicted that the hypothesis does not reflect a chance
correlation.
Fischer’s randomization test
In order to check the strong relationship between the chemical
compound and the biological activity of the training set molecules
the Fischer’s randomization test was performed. In this method, the
biological activity was scrambled and assigned new values. After
that the pharmacophore hypotheses were generated with new
values of activity using the original pharmacophoric features and
constraints used to generate the original pharmacophore
hypotheses. If the randomization run generates improved
correlation coefficient and/or better statistical parameters than the
original hypothesis may be considered to be developed by chance.
Different number of spreadsheets are generated based on the
statistical significance randomization run. The statistical significance
is given by following equation.
[1 (1 )/ ]Significance a b
(1)
Where, a denotes the number of hypotheses with a total cost less
than the best hypothesis, whereas b implies a collection of
HypoGen runs and random runs. For example, at 95% confidence
level total number of random spreadsheet are generated as 19 (b =
20) and each generated spreadsheet is submitted to HypoGen using
the same parameters as the initial run. In the present study, the
developed pharmacophore model was tested at 95% confidence
level which produced 19 spreadsheets.
Test set prediction
The ability to judge the external predictivity of the model beyond
the training set molecules is an important step which verifies ability
of prediction of test compounds. Accordingly, in the present work
410 test compounds were predicted using the developed
pharmacophore model by Ligand Pharmacophore Mapping
protocol in DS, given in Table S1 (Supplementary file). Quality of
prediction of the pharmacophore model was adjudged best on
statistical parameters, R
2
pred
(correlation coefficient) and s
p
(error of
prediction)
39, 40
.
The R
2
pred
value depends on the average experimental activity of the
training set molecules. Since these parameters depend on average
value, it might be achieved for compounds with a wider range of
activity value, but this may not be guaranteed that the estimated
activity values are very close to those experimental activity. As a
result, instead of a good overall correlation being maintained, there
is chance of a significant numerical difference between the two
values. In order to better indicate the predictivity of the
pharmacoinformatics model, modified r
2
[r
2
m(test)
]
41, 42
values were
calculated (threshold value=0.5).
Decoy set
In order to check the efficiency of the screening capacity of the
selected pharmacophore model, it was validated using the decoy
set approach. Decoy set method checks how the model can select
active molecules over inactive molecules on screening with
amalgamated active and inactive molecules. In this purpose a set
decoy was generated by DecoyFinder1.1
43
. Decoys physically
resemble active inhibitors but differ chemically to avoid bias in the
enrichment factor calculation. Based on five parameters decoys
were selected and these included molecular weight, number of
rotational bonds, hydrogen-bond donor count, hydrogen-bond
acceptor count, and the octanol–water partition coefficient of the
active inhibitors. In order to discriminate decoys and active
inhibitors chemically, the MACCS fingerprints were calculated
according to the maximum Tanimoto coefficient values. The
screening dataset consisting of 80 active HIV-integrase inhibitors
and 900 decoys obtained from DecoyFinder was used for queries of
the selected hypothesis. Different statistical parameters including
the accuracy and enrichment factor (EF) were calculated to validate
the pharmacophore model. Screened molecules based on
pharmacophore model were ranked on basis of fit value. For the
assessment of effectiveness of screening the enrichment calculation
of the dataset was performed. The EF implies total known active
inhibitors retrieved from the part of screened database. In the
current study, EF (1%) was calculated from the top 1% hits. Another
parameter Boltzmann-enhanced discrimination of receiver
operating characteristic (BEDROC) which gives the significance of
the dataset screening was also calculated. The BEDROC is a
comprehensive form of receiver operating characteristic (ROC),
which recognises problems in the screening method. Calculation of
the enrichment factor and BEDROC are described by Bhayye et. al
44
.
Virtual screening
5