Arch. Pharm. Res. DOI 10.1007/s12272-013-0313-1

RESEARCH ARTICLE

In silico study on indole derivatives as anti HIV-1 agents: a combined docking, molecular dynamics and 3D-QSAR study Anand Balupuri • Changdev G. Gadhe • Pavithra K. Balasubramanian • Gugan Kothandan Seung Joo Cho

•

Received: 28 September 2013 / Accepted: 3 December 2013 Ó The Pharmaceutical Society of Korea 2013

Abstract The HIV-1 envelope glycoprotein gp120 plays a vital role in the entry of virus into the host cells and is a potential antiviral drug target. Recently, indole derivatives have been reported to inhibit HIV-1 through binding to gp120, and this prevents gp120 and CD4 interaction to inhibit the infectivity of HIV-1. In this work, molecular docking, molecular dynamics (MD) and three-dimensional quantitative structure–activity relationship studies were carried out. Molecular docking studies of the most active and the least active compounds were performed to identify important residues in the binding pocket. We refined the docked poses by MD simulations which resulted in conformational changes. After equilibration, the structure of the ligand and receptor complex was stable. Therefore, we just took the last snapshot as the representative binding pose for this study. This pose for the most active inhibitor was used as a template for receptor-based alignment which was subsequently used for comparative molecular field analysis. Resultant 3D contour maps suggested smaller Anand Balupuri and Changdev G. Gadhe have contributed equally to this work.

Electronic supplementary material The online version of this article (doi:10.1007/s12272-013-0313-1) contains supplementary material, which is available to authorized users. A. Balupuri  C. G. Gadhe  P. K. Balasubramanian  G. Kothandan  S. J. Cho (&) Department of Bio-New Drug Development, College of Medicine, Chosun University, 375 Seosuk-dong, Dong-gu, Gwangju 501-759, Republic of Korea e-mail: [email protected] S. J. Cho Department of Cellular and Molecular Medicine, College of Medicine, Chosun University, Gwangju 501-759, Republic of Korea

substituents are desirable at the 7-position of indole ring to avoid steric interactions with Ser375, Phe382 and Tyr384 residues in the active site. These results can be exploited to develop potential leads and for structure-based drug design of novel HIV-1 inhibitors. Keywords HIV-1 entry inhibitors  gp120  CoMFA  Molecular docking  Molecular dynamics simulation  MM–GBSA

Introduction According to the report of Joint United Nations Programme on HIV/AIDS (http://www.unaids.org/en/media/unaids/con tentassets/documents/epidemiology/2012/gr2012/20121120_ FactSheet_Global_en.pdf), approximately 34.0 million people worldwide were living with HIV at the end of 2011. Also, infection by the HIV virus and the subsequent progression to acquired immunodeficiency syndrome (AIDS) caused death of 1.7 million people worldwide in 2011. This report highlights the urgent need to develop medical protection, prevention and treatments for HIV/AIDS. Therefore, in the past years, intensive efforts have been devoted towards the discovery of novel anti HIV-1 agents. HIV infection is mediated by a series of attachment events. They are initiated by the viral envelope glycoproteins. These glycoproteins are organized into oligomeric, apparently trimeric spikes which are displayed on the surface of the virion (Kowalski et al. 1987; Lu et al. 1995). Every trimeric envelope spike is comprised of three gp120 envelope glycoproteins and three gp41 transmembrane proteins (Blacklow et al. 1995). The HIV-1 entry is a multistep process that includes attachment, CD4 binding, coreceptor binding and membrane fusion (Wyatt and

123

A. Balupuri et al.

Sodroski 1998; Bosch et al. 1989; Brasseur et al. 1988; Helseth et al. 1990; Dorr et al. 2005; Kilby et al. 1998; Teixeira et al. 2009, 2011, 2013; Zhan et al. 2010; Kadow et al. 2006). The binding of viral envelope protein gp120 to the host cell receptor CD4 is the first step in the entry of the virus into target cells. CD4 is a glycoprotein expressed on T-lymphocytes as well as several other cells of the immune system (Berger et al. 1999). This interaction triggers a conformational change in gp120 that facilitate subsequent binding of gp120 to its cellular co-receptors CCR5 or CXCR4. Once the co-receptor binding takes place, gp41 rearranges to form the six-helix bundle which inserts into the host cell membrane. This terminates in the fusion of viral and host membranes by a still unresolved mechanism. Each of these steps as well as the overall process offers new potential targets for drug discovery (Blacklow et al. 1995; Dalgleish et al. 1984; Lasky et al. 1987). Inhibition of the initial entry event of the virus into host cells appears to be a significant, however challenging approach for the prevention or treatment HIV-1 infection and AIDS. Undoubtedly, there have not been extensive studies on developing viral entry inhibitors that target this step (Blacklow et al. 1995). This study is to elucidate the structural features that govern the inhibitory activity of HIV-1 entry inhibitors. In the current work, a set of in silico methods such as molecular docking, MD simulations and three-dimensional quantitative structure–activity relationship (3D-QSAR) studies have been performed.

Molecular alignment for 3D-QSAR analysis Molecular alignment is the most important step of the 3DQSAR model development. The accuracy of comparative molecular field analysis (CoMFA) model prediction and its reliability of contour maps in 3D space strictly depend upon the alignment of the molecules. The most active compound (number 49 in Table 1) in the dataset was selected as a template molecule for sketching and alignment of all other compounds. In order to derive the statistically significant CoMFA models, two different superimposition schemes were investigated. The first scheme was ligand-based superimposition. In this approach, lowest energy conformation of compound 49 derived from systematic conformational search was considered as template for sketching the rest of the molecules. The molecules were constructed by changing the substituents in the template structure. The common substructurebased alignment method was used to align all the inhibitors. The molecular superposition of all data set molecules are shown in Fig. 1a. The second scheme corresponds to the receptor based superimposition which is shown in Fig. 1b. According to this strategy, alignment was generated on the basis of obtained conformation of template molecule derived from docking and 10 nanoseconds (ns) MD simulation. All molecules were re-sketched on the basis of this template structure and minimized in the binding pocket of the receptor. Molecular docking

Materials and methods Data set A data set of 68 indole-based inhibitors along with their biological activities was collected from the literature reported by the same research group (Yeung et al. 2013a, b). The inhibitory activities of HIV-1 inhibitors were expressed as EC50 values. The reported values were converted into pEC50 (-logEC50) for 3D-QSAR study. Chemical structures and their corresponding biological activities are shown in Table 1. Structures of all compounds were sketched using SYBYL 8.1 package (Tripos Inc., St. Louis, MO, USA). Energy minimization of all the molecules was performed using the Tripos molecular mechanics force field (Clark et al. 1989), the Gasteiger– Hu¨ckel charge (Gasteiger and Marsili 1980) with a distance-dependent dielectric and Powell conjugate gradient algorithm (Powell 1977). The minimization process was terminated when the energy gradient convergence criterion of 0.05 kcal mol-1 was attained or when the 10,000 steps minimization cycle was surpassed.

123

Docking studies were performed to identify the probable bioactive conformations of these compounds and to investigate the protein–ligand interactions. For this goal, the most and least active compounds from the data set were docked in the active site of protein. Until recently, no crystal structure of gp120 in complex with indole derivatives has been reported. In this study, a recent crystal ˚ resostructure of HIV-1 gp120 (PDB code: 4DKR, 1.8 A lution) was retrieved from RCSB protein databank. AutoDock 4.2.5.1 suite (Morris et al. 2009) was used as molecular-docking tool, to perform the docking studies. Prior to docking, all the hetero-atoms including co-crystallized ligand as well as water molecules and ions were removed from the protein structure. Polar hydrogen atoms were added and Gasteiger charges were assigned to the protein structure. The resulting target protein structure was used in molecular docking studies. The ligands were prepared for molecular docking by merging non-polar hydrogens, assigning Gasteiger charges, defining the rotatable bonds and assigning AutoDock type to each atom using AutoDock Tools (ADT) 1.5.6 (Sanner 1999).

In silico study on indole derivatives as anti HIV-1 agents Table 1 Data set of indole derivatives employed in the 3D QSAR study

Core structure

Compound 2 Substituent

Activity

Compound R

X

Y

Z

EC50 (nM)

pEC50

-

F

H

C

1.240

8.907

-

-

-

-

0.750

9.125

F

H

C

2.030

8.693

F

H

C

0.520

9.284

F

H

C

407.000

6.390

F

H

C

2.680

8.572

F

H

C

3.830

8.417

F

H

C

2.690

8.570

F

H

C

18.820

7.725

F

H

C

12.830

7.892

F

H

C

315.700

6.501

F

H

C

7.950

8.100

F

H

C

11.240

7.949

F

H

C

109.100

6.962

F

H

C

6.730

8.172

1 2 H2N

3

O HN

4 O

5a,b

N O

O

6

N

N Ph

O

MeO

7

H N O

H2N

8 9

H N O H N

MeO

O

10

H N

H2N

O H N

N

11

O H N

N

12 O

O

O

13

N

HN O

HN

14

15

O

H N

HN

O

N

123

A. Balupuri et al. Table 1 continued

Substituent

Activity

Compound R

16

O

17

S

HN

HN

Y

Z

18

O

F

H

C

1.300

8.886

O

F

H

C

1.060

8.975

F

H

C

1.080

8.967

F

H

C

2.640

8.578

F

H

C

0.290

9.538

OMe

H

C

0.760

9.119

OMe

H

C

0.760

9.119

F

Me

C

0.140

9.854

F

Me

C

0.750

9.125

F

H

N

7.590

8.120

F

Me

N

1.420

8.848

F

H

C

3.370

8.472

F

H

C

9.090

8.041

F

H

C

0.650

9.187

F

H

C

0.400

9.398

F

H

C

0.190

9.721

F

H

C

3.070

8.513

O O S HN O

19

pEC50

O

O Ph S HN O

O

H2N

21

EC50 (nM)

H N

H N

20

X

O HN

22 O H2N

23

O

HN

24 O H2N

25

O H2N

26

O

27

N

HN O

28

HN N O HN

29

N O

30

H N N N N

HN O

S

31

N

HN N O

32

N

HN N O

123

In silico study on indole derivatives as anti HIV-1 agents Table 1 continued

Substituent

Activity

Compound R

X

Y

Z

EC50 (nM)

pEC50

F

H

C

0.940

9.027

F

H

C

0.020

10.699

F

H

C

0.057

10.244

F

H

C

0.006

11.237

F

H

C

0.320

9.495

F

H

C

0.040

10.398

F

H

C

0.210

9.678

H

H

C

153.000

6.815

F

H

C

0.430

9.367

F

H

C

0.040

10.398

F

H

C

0.120

9.921

F

H

C

0.060

10.222

F

H

C

0.060

10.222

F

H

C

0.010

11.000

F

H

C

0.007

11.174

F

H

C

0.050

10.301

F

H

C

0.005

11.268

F

H

C

0.060

10.222

F

H

C

0.050

10.301

F

H

C

0.050

10.301

N O HN

33

O

S

34

a

HN

N

O S

35

HN

N

O

S

36

HN

N

O

S

37

HN

N

O NH

38

HN

N

O

S

39

HN

N

Ph

O

40a

-

41 42

MeO

43 44

N N

45 46

S

S

N

47

S N

48 49b 50 51 52

O O N N

N H N

O N N N O

123

A. Balupuri et al. Table 1 continued

Substituent

Activity

Compound R

X

Y

Z

EC50 (nM)

pEC50

F

H

C

0.030

10.523

F

H

C

0.080

10.097

F

H

C

0.880

9.056

F

H

C

0.530

9.276

F

H

C

0.060

10.222

F

H

C

0.094

10.027

F

H

C

3.550

8.450

F

H

C

0.008

11.086

F

H

C

0.040

10.398

F

H

C

0.300

9.523

F

H

C

0.100

10.000

F

H

C

0.100

10.000

F

H

C

0.160

9.796

OMe

H

C

0.160

9.796

F

Me

C

0.040

10.398

F

Me

C

0.070

10.155

N N

53

N H N N

54

HN N N

55

O N

56

N H N

57 58

O N N N O MeO N

59

O N

H2N

60

N O N

HN

61

N O N

N

62

N O N

HN

63

O

N O N

NH2

64

N

O

O N N N

65 66 67 68

N H N

O N N O N N N

HN N

a

Compounds are considered as outliers

b

Compounds 5 and 49 are referred as Low and High, respectively, in the manuscript

˚ grid points with A 3D grid box with 40 9 40 9 40 A ˚ grid spacing of 0.375 A was created by the AutoGrid algorithm. Center of the grid was set using Ca atom of the active site residue Trp427 of the macromolecule.

123

Lamarckian genetic algorithm (LGA) method (Morris et al. 1998) was employed to search for the globally optimized conformations. With remaining parameters set as default values, each docking experiment was derived from 100

In silico study on indole derivatives as anti HIV-1 agents

Fig. 1 Alignment of data set molecules based on common substructure using compound 49 as a template. a Ligand-based alignment. b Receptor-based alignment. (Color figure on the web)

different runs that were set to terminate after a maximum of 250,000 energy evaluations and eventually generated 100 docked conformations for each ligand. Docking results of both experiments were analyzed by the cluster analysis. For both ligands, binding modes with lower binding free energies and larger cluster populations than the rest modes were selected for further analysis.

equilibrate the solvent density. Finally, a production run for 10 ns was performed at 300 K and 1.0 atm pressure with a time step of 2 fs. Trajectory analysis and the binding free energy calculations were carried out using the PTRAJ and MM–GBSA modules in AMBER package.

Molecular dynamics (MD) simulation

The binding free energies were estimated by MM–GBSA method implemented in AMBER 12. This procedure employs an all-atom method that uses molecular mechanics with continuum solvent and the OPLS_2005 force field to compute the binding free energies. A total of 400 snapshots were extracted from the last 2 ns of MD trajectories and used for the energy calculations. In MM–GBSA, binding free energy (DGbinding) is computed for each molecular species (complex, receptor and ligand) as the difference between the free energy of the complex (Gcomplex) and the free energy of the receptor (Greceptor) plus the ligand (Gligand).
i:e:; DGbinding ¼ Gcomplex  Greceptor þ Gligand

MD simulations were carried out to refine the docking results using AMBER 12 package (Case et al. 2012). The two docked complexes of gp120 with the most and least active compounds were used as the initial structures for MD simulations. The FF99SB AMBER force field (Hornak et al. 2006) was used for the protein while ligand parameters were generated by the general amber force field (GAFF) (Wang et al. 2004). Ligand–protein complex was immersed in a global box of TIP3P water molecules with a ˚ . One sodium counter ion was added margin distance of 8 A to the complex to ensure overall charge neutrality of the system. The particle mesh Ewald (PME) (Essmann et al. 1995) was applied to deal with long-range electrostatic ˚ cut-off distance was used for the interactions and 12 A non-bonded interactions. With a time step of 2.0 femtoseconds (fs), the SHAKE algorithm was adopted to constrain the bond lengths. Periodic boundary conditions were employed to avoid edge effects in all calculations. Subsequently, system was energetically minimized by the steepest-descent and conjugate gradient methods to eliminate possible bad contacts. The two stage equilibration was performed on the system. Firstly, the system was gradually heated from 0 to 300 K, over 50 picoseconds (ps) of simulation time. Further, the system was equilibrated for 100 ps using pressure and temperature control in order to

Binding free energy calculations

Each term can be computed as follows: DG ¼ DGMM þ DGsol TDS DGMM ¼ DGele þ DGvdw DGsol ¼ DGele;sol þ DGnonpol;sol where DGbinding is binding free energy. TDS represents entropy terms. The molecular mechanics energy composes of electrostatic interactions (DGele) and van der Waals interactions (DGvdw). Solvation free energy includes polar (DGele,sol) and the nonpolar contributions (DGnonpol,sol). MM–GBSA determines binding free energy in the similar manner as molecular mechanics–Poisson Boltzmann

123

A. Balupuri et al.

surface area (MM–PBSA). However, MM–PBSA uses the Poisson-Boltzmann model while MM–GBSA uses generalized Born model to calculate the electrostatic solvation energy. A recent report revealed that in most systems, MM/ GBSA performed better than MM/PBSA in predicting relative binding free energies (Hou et al. 2011). Therefore, MM–GBSA calculations were performed to identify the most favorable binding mode. Calculation of CoMFA descriptors The steric and electrostatic CoMFA potential fields for all molecules were calculated at each lattice intersection of a ˚ . The grid box dimensions regularly spaced grid of 2.0 A were determined automatically in such a way that region ˚ in each direction boundaries were extended beyond 4.0 A from coordinates of each molecule. The CoMFA fields were calculated separately for each molecule using a default probe atom, a sp3 hybridized carbon atom with ?1 ˚ . The Lencharge and a van der Waals radius of 1.52 A nard–Jones potential terms and Coulombic terms signifying steric and electrostatic fields respectively, were calculated using Tripos force field. Energy cut-off value of 30.0 kcal mol-1 was selected for both steric and electrostatic fields. A partial least squares (PLS) regression (Wold et al. 1984) was implemented to obtain a correlation between the structural parameters (CoMFA descriptors) and experimental biological activities (pIC50 values). The cross-validation analysis was performed using the leave-

one-out (LOO) method to check the predictivity of the derived model and to determine the optimum number of components (ONC) with minimum standard error of prediction (SDEP). In LOO method, one molecule is omitted from the data set and its activity is predicted by the model derived from the rest of the data set molecules. The minimum column filtering was set to 2.0 kcal mol-1 during model generation to accelerate the regression analysis and to improve the signal-to-noise ratio by omitting those lattice points whose energy variation was below this threshold. The cross-validated q2 with ONC and the lowest SDEP was considered for further analyses. Finally, non-cross-validated analysis was performed to calculate conventional Pearson correlation coefficient (r2), standard error of estimate (SEE) and Fischer’s statistical value (Fratio) using the previously obtained ONC. The region focusing method was applied to refine the conventional CoMFA models.

Results and discussion Molecular docking In order to illustrate the binding pose of HIV-1 entry inhibitors to gp120 and to predict binding interactions, the most potent compound 49 (High) and least active compound 5 (Low) were selected as representatives. Both compounds were docked in the binding pocket of gp120. Docked conformations of the compounds are shown in

Fig. 2 Docked conformations of a High and b Low in the active site of gp120. Active site residues are shown as lines and ligand is displayed as stick model. Non-polar hydrogens are omitted for clarity. (Color figure on the web)

123

In silico study on indole derivatives as anti HIV-1 agents

Fig. 2a, b. Both compounds were well positioned in the binding pocket in the same way and displayed a similar binding mode. Docking results demonstrated that indole ring was placed deep into the gp120 cavity enclosed by hydrophobic residues whereas the piperazine ring was docked at the entrance of binding pocket lined by Ile371, Ile424, Asn425, Met426 and Trp427 residues. The hydrophobic site was composed of the Val255, Phe376, Phe382, Tyr384, Ile424, Met426, Trp427, Gly473 and Ile475 residues. This observation was in accordance with the previous studies (Xiang et al. 2002; Gadhe et al. 2012; Teixeira et al. 2013). As shown in Fig. 2a, High was inserted deep into crevices of target protein with oxazole ring directed towards the Val255, Ser256 and Asn377. Hydrogen bond interactions were missing for the ligand in the binding site. High showed both hydrophobic and hydrophilic interactions in the binding site. It was noticed that the residues Ser256, Ser375 and Asn377 were involved in hydrophilic interaction whereas Phe376 and Phe382 were participating in hydrophobic interactions. As can be seen in Fig. 2b, hydrogen bond interactions were absent for Low. This compound was accommodated in the active site primarily through interactions of tertiary amide group with Trp112, Val255, Ser375 and Asn377 residues. Docking simulations contributed to the identification of the some key residues in the active site such as Trp112, Val255, Ser256, Ser375, Phe376, Asn377, Phe382, Asn425, Met426 and Trp427. These residues were the main contributors to the interactions between ligand and receptor. MD simulations Although docking analysis can provide a good starting point for predicting binding pose but solvent effect and ligand-induced conformational changes are not fully taken into account. Therefore, 10 ns MD simulations were performed on both the docked complexes to investigate conformational changes and to further explore ligand–receptor interactions. The dynamic stability of the MD trajectories was monitored by the root mean square deviations (RMSD) values of backbone atom with respect to the starting structure. Temperature, pressure and potential energy of the systems were also monitored during the simulations and curves are depicted in Fig. S1a–c, respectively (supplementary data). It can be noted that both complexes share a similar RMSD behavior. The graph in Fig. 3a shows that both the systems attained equilibrium within 5,000 ps and remained stable throughout simulations. The mean RMSD values of High-gp120 and Low-gp120 systems were 1.70 ˚ , respectively, and the relative RMSD fluctuaand 1.72 A tions were very small. Figure 3b displays RMSD of two

inhibitors which demonstrates that Low was less stable as compared to the High during the simulations. It can be seen that average RSMD fluctuations for High and Low ˚ , respectively. For both systems, it were 0.83 and 1.52 A was observed that MD simulated and docked structures of complexes were located in the same binding pocket. Some positional deviations were observed between the docked and MD simulated structures for Low-gp120 system. One can notice that orientation of the Low was changed at the entrance of binding cavity after MD simulation. The piperazine and phenyl ring of this compound shifted to a direction opposite to that of docked orientation. The orientation as well as conformation of High remained almost same after the MD simulation. The binding pose obtained for this compound was in agreement with reported study (Kong et al. 2006). Furthermore, analysis of root-meansquare fluctuation (RMSF) versus the residue number for complexes is shown in Fig. 3c. Higher fluctuations were observed for Low-gp120 complex as compared to Highgp120. These results suggest that binding pose of Low was not stable. Thus, MD simulated pose of High was selected for the receptor-based CoMFA. The hydrogen bond interaction is one of the most essential forces that maintain the binding between ligand and receptor. For High-gp120 system, two hydrogen bond interactions were observed as shown in Fig. 3d. The NH group of the indole ring established a hydrogen bond with side chain OH group of Ser375 and piperazinamide oxygen formed another hydrogen bond with side chain amine group of Asn478. In Low-gp120 system, benzoyl group formed a hydrogen bond with amine group of Asn425 located at the entrance site of the binding pocket as displayed in Fig. 3e. Binding free energy analysis The computed binding free energies (DGbinding) of both systems are listed in Table 2. The High-gp120 system was found to be energetically (DGbinding) more stable than Lowgp120 system. Data shows that van der Waals interactions made major contribution to the binding and suggested that systems were having good hydrophobic contacts. The nonpolar solvation contribution terms which refers to the burial of solvent accessible surface area upon binding, provided slightly favorable contribution to the binding. The unfavorable polar solvation contribution terms opposed interactions within the protein–ligand complexes. The energy component analysis displayed that sum of the van der Waal interaction energy (-50.280 kcal mol-1) and electrostatic energy (-14.927 kcal mol-1) in High-gp120 was higher than the Low-gp120 (-49.443 and -8.763 kcal mol-1). Also, the cost of solvation free energy (DGnonpol,sol ? DGele,sol) of Low-gp120 was high (28.365 kcal mol-1) as

123

A. Balupuri et al.

123

In silico study on indole derivatives as anti HIV-1 agents b Fig. 3 a Backbone-RMSD, b inhibitors-RMSD and c RMSF curves

for both the simulated systems. MD-simulated structures of d High and e Low with the binding site. (Color figure on the web)

Table 4 Experimental and predicted pEC50 values of data set generated by CoMFA models along with their residuals Compounds

pEC50

Table 2 Binding free energy (kcal mol-1) of both systems Energy

High-gp120

Ligand-based CoMFA

Receptor-based CoMFA

Predicted

Residual

Predicted

Residual 0.036

Low-gp120

1

8.907

8.811

0.096

8.871

-49.4

2

9.125

8.878

0.247

8.826

0.299

-14.9

-8.8

3

8.693

9.482

-0.790

9.179

-0.487

-5.9

-5.5

4

9.284

8.970

0.314

9.090

0.194

30.2

33.9

6

8.572

8.764

-0.192

8.628

-0.056

-65.2

-58.2

7

8.417

8.571

-0.154

9.037

-0.620

28.4

8

8.570

8.873

-0.303

9.764

-1.194

-25.8

-20.3

9

7.725

8.797

-1.072

8.206

-0.481

-15.1

-9.5

10

7.892

7.393

0.499

8.272

-0.380

11

6.501

7.485

-0.984

7.107

-0.606

12

8.100

8.112

-0.012

7.750

0.350

13

7.949

8.395

-0.446

7.745

0.204

14

6.962

7.405

-0.443

7.496

-0.534

15 16

8.172 8.886

8.064 8.395

0.108 0.491

8.337 8.586

-0.165 0.300

17

8.975

8.524

0.451

8.745

0.230

18

8.967

8.691

0.276

8.650

0.317

19

8.578

8.332

0.246

8.257

0.321

20

9.538

9.494

0.044

9.046

0.492

Table 3 Statistical results of CoMFA models

21

9.119

9.276

-0.157

9.088

0.031

Parameters

22

9.119

9.112

0.007

9.025

0.094

23

9.854

9.449

0.405

9.207

0.647

24

9.125

9.286

-0.161

9.505

-0.380

25

8.120

9.077

0.025

8.671

-0.551

26

8.848

8.730

0.118

8.656

0.192

27

8.472

8.252

0.220

8.232

0.240

28

8.041

8.488

-0.447

8.092

-0.051

29

9.187

9.524

-0.337

9.850

-0.663

30 31

9.398 9.721

9.518 9.220

-0.120 0.501

9.984 9.862

-0.586 -0.141

DGavdw DGbele DGcnonpol,sol DGdele,sol DGeMM DGfsol g

-50.3

TDS

DGhbinding a b c

24.3

van der Waals interaction energy Electrostatic interaction energy Nonpolar contribution

d

Polar contribution

e

Molecular mechanics energy

f

Solvation free energy

g

Entropy

h

Estimated binding free energy

Ligand-based CoMFA

Receptor-guided CoMFA

Without region focusing

With region focusing

Without region focusing

With region focusing

q2a

0.544

0.587

0.486

0.506

2b

0.844

0.815

0.815

0.802

SEPc

0.711

0.670

0.751

0.736

d

0.415

0.449

0.451

0.466

r

SEE

Feratio ONCf

52.385 6

52.03 5

65.967 4

60.63 4

32

8.513

8.218

0.295

8.621

-0.108

a

Cross validated correlation coefficient

33

9.027

9.428

-0.401

10.144

-1.117

b

Non-cross validated correlation coefficient

35

10.244

10.007

0.237

9.930

0.314

c

Standard error of prediction

36

11.237

10.138

1.099

10.499

0.738

d

Standard error of estimation Fischer’s statistical value

37

9.495

10.206

-0.711

9.655

-0.160

38

10.398

10.344

0.054

9.709

0.689

Optimum number of components

39

9.678

9.333

0.345

9.827

-0.149 -0.381

e f

-1

compared to the High-gp120 (24.321 kcal mol ). Thus, the binding of the High-gp120 was stronger than that of the Low-gp120 complex. 3D-QSAR models CoMFA models were developed for indole derivatives as HIV-1 entry inhibitors using ligand-based and receptor-

41

9.367

9.611

-0.245

9.747

42

10.398

9.841

0.557

9.769

0.629

43

9.921

10.165

-0.244

9.977

-0.056

44

10.222

10.530

-0.308

9.904

0.318

45

10.222

10.019

0.203

10.491

-0.269

46

11.000

11.083

-0.083

10.429

0.571

47

11.174

10.807

0.367

10.698

0.476

48 49

10.301 11.268

10.670 10.899

-0.369 0.369

10.874 10.624

-0.573 0.644

123

A. Balupuri et al. Table 4 continued Compounds

pEC50

Ligand-based CoMFA

Receptor-based CoMFA

Predicted

Predicted

Residual

Residual

50

10.222

10.004

0.218

9.999

0.223

51

10.301

10.329

-0.028

10.442

-0.141

52

10.301

10.337

-0.036

9.790

0.511

53 54

10.523 10.097

10.202 10.528

0.321 -0.431

10.559 10.591

-0.036 -0.494

55

9.056

9.761

-0.706

9.535

-0.480

56

9.276

9.203

0.073

9.523

-0.247

57

10.222

10.395

-0.173

10.306

-0.084

58

10.027

10.123

-0.096

9.268

0.759

59

8.450

8.613

-0.163

7.968

0.482

60

11.086

10.377

0.709

10.802

0.284

61

10.398

9.445

0.953

10.512

-0.114

62

9.523

8.790

0.733

9.270

0.253

63

10.000

10.351

-0.351

9.364

0.636

64

10.000

9.744

0.256

9.604

0.396

65

9.796

9.516

0.280

9.825

-0.029

66

9.796

10.316

-0.520

10.456

-0.660

67

10.398

10.646

-0.248

10.394

0.004

68

10.155

10.539

-0.384

10.034

0.121

based alignment schemes. The developed models were analyzed by several statistical parameters which are summarized in Table 3. The data set was not divided into training set and test set during model generation. From the data set, three compounds (5, 34 and 40) were removed as outliers. We inspected these molecules and found that compound 40 was structurally different from the rest of the molecules in the dataset. This compound lacked substitutions at 4-position (X) as well as 7-position (R) of the indole ring. Compound 5 was found to be least active in the series with EC50 of 407 nM and very different from others in terms of activity. Models generated without compounds 5 and 40 using both ligand-based and receptor-based

Fig. 4 Scatter plot of the experimental versus predicted activities. a Ligand-based CoMFA and b Receptor-based CoMFA. (Color figure on the web)

123

alignment schemes predicted compound 34 with a residual error larger than 1 logarithm unit. Thus, compound 34 was also considered as outlier. So, the final models were generated on the basis of 65 molecules. Removal of outliers produced statistically reasonable results for both ligandbased (q2 = 0.544, r2 = 0.844) and receptor-based (q2 = 0.486, r2 = 0.815) models. Subsequently, these models were optimized through the region focusing which further improved statistical parameters. This approach produced more reliable CoMFA region-focusing models for ligand-based (q2 = 0.587, r2 = 0.815) and receptorbased (q2 = 0.506, r2 = 0.802) alignment schemes. It can be seen in Table 3, that application of region focusing improved the internal predictability (q2) from 0.544 to 0.587 of the ligand-based model and 0.486 to 0.506 of the receptor-based model but negligibly reduced their r2 values. Although, the final region-focused model resulted from the ligand-based alignment has a slightly better statistics, receptor-guided alignment provides information about the inhibitors and receptor interactions. Thus, we focused on CoMFA region-focusing model obtained using receptor-based alignment for further discussion. Contour maps of the model were taken as representative to explain factors affecting activities of inhibitors. The experimental and predicted pEC50 values along with residual values for data set are presented in Table 4 and the relationship between experimental and predicted pEC50 values is depicted in Fig. 4a, b. The contour plots of the CoMFA steric and electrostatic fields are displayed in Fig. 5a, b, respectively. To aid in visualization, the most potent High is overlaid with the maps. In steric contour map, green contour depicts the sterically favorable regions and yellow contour represents unfavorable region. From Fig. 5a, it can be seen that the green contour map was observed in the vicinity of oxazole ring. It indicates that the smaller bulky substituents are favorable for improving inhibitory potency. In case of compounds 3 (pEC50 = 8.693) and 4 (pEC50 = 9.284) smaller substituents like amide and

In silico study on indole derivatives as anti HIV-1 agents

Fig. 5 CoMFA contour plots with the combination of High as a reference. a Steric contour map. Green contours indicate regions where bulky groups increase activity and yellow contours indicate

regions where bulky groups decrease activity. b Electrostatic contour map. Blue contours indicate regions where positive charges increase activity. (Color figure on the web)

secondary methyl amide at 7-position of indole ring are well tolerated and found to retained the activity. However, the tertiary N containing compounds 5 (pEC50 = 6.39) and 11 (pEC50 = 6.50) were found to be less tolerated in the binding site. It was observed that the binding pocket around oxazole ring is sterically restricted by Ser375, Phe382 and Tyr384 residues and thus, smaller substituents are favorable for improving inhibitory potency. However, aromatic rings with the heteroatom (compounds 15–18, pEC50 = 8.172–8.975) are well tolerated in the binding site than just aromatic ring (compound 14, pEC50 = 6.962). The reason for this affinity improvement could be the formation of hydrogen bonds between the heteroatom containing aryl ring and active site residues which may not be possible in case of the aromatic ring. Two small yellow contour maps were observed near the Ser256 and Tyr384, which indicates that the substitution penetrating that direction are not favorable for higher inhibitory potency. In case of compounds 9–13 (pEC50 = 6.501–8.1), more extended amide substituents could sterically interacts with pocket wall residues (Ser256, Thr257 and Tyr384) and appeared to be medium active. In electrostatic contour map, the blue and red contour represents regions where electropositive and electronegative groups are favorable for inhibitory potency. From Fig. 5b, it can be seen that two medium sized blue contour maps (electropositive favorable) were observed in the vicinity of oxazole ring of template molecule. Comparison between compound 41 (pEC50 = 9.367) and compound 48 (pEC50 = 11.27) suggests that the aromatic substituents at 7-position of indole ring is less favorable than the heteroaryl group. Electropositive blue contour suggest that the polar

group are favorable at this position for improved inhibitory potency, because surrounding pocket is lined by the polar residues. It can be also suggested that the interaction between heteroaryl moiety and pocket residues are highly specific and orientation driven. This is the reason why compounds 45 (pEC50 = 10.22) and 46 (pEC50 = 11) had difference in the activity. In case of 1,2,4-oxadiazole substituted derivatives 60 (pEC50 = 11.09), 61 (pEC50 = 10.4) and 62 (pEC50 = 9.52), it can be seen that the activity is reduced from the primary \ secondary \ tertiary amine. The hydrophilic nature of primary amine group is more favorable than tertiary amine to improve the potency, because the pocket is surrounded by the hydrophilic residues. From the obtained QSAR results, it can be concluded that the binding pocket is sterically restricted and the smaller substituents with balanced steric and electrostatic property are highly desirable at 7-position of indole ring to improve the inhibitory potency. These observations are in agreement with the previous study (Teixeira et al. 2009).

Conclusions In the present work, molecular docking, MD simulations and 3D-QSAR studies were carried out. These approaches were applied to elucidate the structural features that govern the activity of HIV-1 entry inhibitors. In addition, this study can also contribute to the better understanding of binding pose between these inhibitors and gp120 glycoprotein. Docking and MD simulations were performed to identify the probable bioactive conformations of these compounds and to investigate the protein–ligand interactions. These simulations

123

A. Balupuri et al.

contributed to the identification of the some important residues in the active site such as Trp112, Val255, Ser256, Thr257, Ser375, Phe376, Asn377, Phe382, Asn425, Met426 and Trp427. These residues play a significant role in the binding of indole based inhibitors to gp120. On the basis of the receptor-based alignment, a reliable region-focusing CoMFA model was obtained. Contour map analysis suggested that the binding pocket is sterically restricted and the smaller substituents with balanced steric and electrostatic properties are highly desirable at 7-position of indole ring to improve the inhibitory potency. All these applied techniques could provide structural information affecting activity of these inhibitors and improve the understanding of ligand– receptor interactions. These results may provide some useful and rational suggestions for further design of novel inhibitors for HIV/AIDS therapy.

Conflict of interest The authors declare that they have no conflicts of interest.

References Berger, E.A., P.M. Murphy, and J.M. Farber. 1999. Chemokine receptors as HIV-1 coreceptors: Roles in viral entry, tropism, and disease. Annual Review of Immunology 17: 657–700. Blacklow, S.C., M. Lu, and P.S. Kim. 1995. Trimeric subdomain of the simian immunodeficiency virus envelope glycoprotein. Biochemistry 34: 14955–14962. Bosch, M.L., P.L. Earl, K. Fargnoli, S. Picciafuoco, F. Giombini, F. Wong-Staal, and G. Franchini. 1989. Identification of the fusion peptide of primate immunodeficiency viruses. Science 244: 694–697. Brasseur, R., B. Cornet, A. Burny, M. Vandenbranden, and J.M. Ruysschaert. 1988. Mode of insertion into a lipid membrane of the N-terminal HIV gp41 peptide segment. AIDS Research and Human Retroviruses 4: 83–90. Clark, M., R.D. Cramer, and V.N. Opdenbosch. 1989. Validation of the general purpose Tripos 5.2 force field. Journal of Computational Chemistry 10: 982–1012. Case, D.A., T.A. Darden, T.E. Cheatham III, C.L. Simmerling, J. Wang, R.E. Duke, R. Luo, R.C. Walker, W. Zhang, and K.M. Merz. 2012. AMBER 12. San Francisco: University of California. Dalgleish, A.G., P.C.L. Beverley, P.R. Clapham, D.H. Crawford, M.F. Greaves, and R.A. Weiss. 1984. The CD4 (T4) antigen is an essential component of the receptor for the AIDS retrovirus. Nature 312: 763–767. Dorr, P., M. Westby, S. Dobbs, P. Griffin, B. Irvine, M. Macartney, J. Mori, G. Rickett, C. Smith-Burchnell, and C. Napier. 2005. Maraviroc (UK-427,857), a potent, orally bioavailable, and selective small-molecule inhibitor of chemokine receptor CCR5 with broad-spectrum anti-human immunodeficiency virus type 1 activity. Antimicrobial Agents and Chemotherapy 49: 4721–4732. Essmann, U., L. Perera, M.L. Berkowitz, T. Darden, H. Lee, and L.G. Pedersen. 1995. A smooth particle mesh Ewald method. Journal of Chemical Physics 103: 8577–8593. Gadhe, C.G., G. Kothandan, T. Madhavan, and S.J. Cho. 2012. Molecular modeling study of HIV-1 gp120 attachment inhibitors. Medicinal Chemistry Research 21: 1892–1904.

123

Gasteiger, J., and M. Marsili. 1980. Iterative partial equalization of orbital electronegativity—a rapid access to atomic charges. Tetrahedron 36: 3219–3228. Helseth, E., U. Olshevsky, D. Gabuzda, B. Ardman, W. Haseltine, and J. Sodroski. 1990. Changes in the transmembrane region of the human immunodeficiency virus type 1 gp41 envelope glycoprotein affect membrane fusion. Journal of Virology 64: 6314–6318. Hornak, V., R. Abel, A. Okur, B. Strockbine, A. Roitberg, and C. Simmerling. 2006. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins: Structure Function, and Bioinformatics 65: 712–725. Hou, T., J. Wang, Y. Li, and W. Wang. 2011. Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations. Journal of Chemical Information and Modeling 51: 69–82. Kadow, J., H.G. Wang, and P.F. Lin. 2006. Small-molecule HIV-1 gp120 inhibitors to prevent HIV-1 entry: An emerging opportunity for drug development. Current Opinion in Investigational Drugs 7: 721–726. Kilby, J.M., S. Hopkins, T.M. Venetta, B. Dimassimo, G.A. Cloud, J.Y. Lee, L. Alldredge, E. Hunter, D. Lambert, and D. Bolognesi. 1998. Potent suppression of HIV-1 replication in humans by T-20, a peptide inhibitor of gp41-mediated virus entry. Nature Medicine 4: 1302–1307. Kong, R., J.J. Tan, X.H. Ma, W.Z. Chen, and C.X. Wang. 2006. Prediction of the binding mode between BMS-378806 and HIV1 gp120 by docking and molecular dynamics simulation. Biochimica et Biophysica Acta (BBA)-Proteins and Proteomics 1764: 766–772. Kowalski, M., J. Potz, L. Basiripour, T. Dorfman, W.C. Goh, E. Terwilliger, A. Dayton, C. Rosen, W. Haseltine, and J. Sodroski. 1987. Functional regions of the envelope glycoprotein of human immunodeficiency virus type 1. Science 237: 1351–1355. Lasky, L.A., G. Nakamura, D.H. Smith, C. Fennie, C. Shimasaki, E. Patzer, P. Berman, T. Gregory, and D.J. Capon. 1987. Delineation of a region of the human immunodeficiency virus type 1 gp120 glycoprotein critical for interaction with the CD4 receptor. Cell 50: 975–985. Lu, M., S.C. Blacklow, and P.S. Kim. 1995. A trimeric structural domain of the HIV-1 transmembrane glycoprotein. Nature Structural & Molecular Biology 2: 1075–1082. Morris, G.M., D.S. Goodsell, R.S. Halliday, R. Huey, W.E. Hart, R.K. Belew, and A.J. Olson. 1998. Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. Journal of Computational Chemistry 19: 1639–1662. Morris, G.M., R. Huey, W. Lindstrom, M.F. Sanner, R.K. Belew, D.S. Goodsell, and A.J. Olson. 2009. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. Journal of Computational Chemistry 30: 2785–2791. Powell, M.J.D. 1977. Restart procedures for the conjugate gradient method. Mathematical Programming 12: 241–254. Sanner, M.F. 1999. Python: A programming language for software integration and development. Journal of Molecular Graphics and Modelling 17: 57–61. Teixeira, C., J.R.B. Gomes, P. Gomes, F.F. Maurel, and F. Barbault. 2011. Viral surface glycoproteins, gp120 and gp41, as potential drug targets against HIV-1: Brief overview one quarter of a century past the approval of zidovudine, the first anti-retroviral drug. European Journal of Medicinal Chemistry 46: 979–992. Teixeira, C., N. Serradji, S. Amroune, K. Storck, C. Rogez-Kreuz, P. Clayette, F. Barbault, and F. Maurel. 2013. Is the conformational flexibility of piperazine derivatives important to inhibit HIV-1 replication? Journal of Molecular Graphics and Modelling 44: 91–103.

In silico study on indole derivatives as anti HIV-1 agents Teixeira, C., N. Serradji, F. Maurel, and F. Barbault. 2009. Docking and 3D-QSAR studies of BMS-806 analogs as HIV-1 gp120 entry inhibitors. European Journal of Medicinal Chemistry 44: 3524–3532. Wang, J., R.M. Wolf, J.W. Caldwell, P.A. Kollman, and D.A. Case. 2004. Development and testing of a general amber force field. Journal of Computational Chemistry 25: 1157–1174. Wold, S., A. Ruhe, H. Wold, and I.W. Dunn. 1984. The collinearity problem in linear regression. The partial least squares (PLS) approach to generalized inverses. SIAM Journal on Scientific and Statistical Computing 5: 735–743. Wyatt, R., and J. Sodroski. 1998. The HIV-1 envelope glycoproteins: Fusogens, antigens, and immunogens. Science 280: 1884–1888. Xiang, S.H., P.D. Kwong, R. Gupta, C.D. Rizzuto, D.J. Casper, R. Wyatt, L. Wang, W.A. Hendrickson, M.L. Doyle, and J. Sodroski. 2002. Mutagenic stabilization and/or disruption of a CD4-bound state reveals distinct conformations of the human

immunodeficiency virus type 1 gp120 envelope glycoprotein. Journal of Virology 76: 9888–9899. Yeung, K.S., Z. Qiu, Q. Xue, H. Fang, Z. Yang, L. Zadjura, C.J. D’Arienzo, B.J. Eggers, K. Riccardi, and P.Y. Shi. 2013a. Inhibitors of HIV-1 attachment. Part 7: Indole-7-carboxamides as potent and orally bioavailable antiviral agents. Bioorganic & Medicinal Chemistry Letters 23: 198–202. Yeung, K.S., Z. Qiu, Z. Yin, A. Trehan, H. Fang, B. Pearce, Z. Yang, L. Zadjura, C.J. D’Arienzo, and K. Riccardi. 2013b. Inhibitors of HIV-1 attachment. Part 8: The effect of C7-heteroaryl substitution on the potency, and in vitro and in vivo profiles of indolebased inhibitors. Bioorganic & Medicinal Chemistry Letters 23: 203–208. Zhan, P., W. Li, H. Chen, and X. Liu. 2010. Targeting protein–protein interactions: A promising avenue of anti-HIV drug discovery. Current Medicinal Chemistry 17: 3393–3409.

123