Traditional Chinese medicine application in HIV: An in silico study
Hung-Jin Huanga†, Yi-Ru Jianb†, Calvin Yu-Chian Chen b, c d, e, f*
a Department of Chinese Pharmaceutical Sciences and Chinese Medicine Resources, China Medical University, Taichung, 40402, Taiwan.
b Department of Biomedical Informatics, Asia University, Taichung, 41354, Taiwan.
c Department of Medical Research, China Medical University Hospital, Taichung, 40447, Taiwan.
d Laboratory of Computational and Systems Biology, China Medical University, Taichung, 40402, Taiwan.
e Department of Biotechnology, Asia University, Taichung, 41354, Taiwan.
f China Medical University Beigang Hospital, Yunlin, 65152, Taiwan.
† equal contribution
*Corresponding author (C. Y. C. Chen).
Tel.: +886-4-2205-2121 ext. 4335
E-mail address: ycc929@MIT.EDU (C.Y.C. Chen) 1
2 3
4 5 6 7 8 9 10 11 12 13 14 15 16 17
Abstract
Viral infection by human immunodeficiency virus (HIV) requires integration of viral DNA with host DNA which involves binding of HIV integrase (IN) with its co-factor len epithelium derived growth factor (LEDGF/p75). Since disrupted binding of IN with LEDGF/p75 inhibits proliferation of HIV, inhibition or denaturation of IN is a possible method for inhibiting HIV replication. D77 is a known drug with demonstrated inhibition against HIV by binding to IN. Herein, we utilized D77 as a control to screen for traditional Chinese medicine compounds that exhibit similar atomic level characteristics. 9-Hydroxy- (10E)-octadecenoic acid and Beauveriolide I were found to have higher Dock Scores to IN than D77 through virtual screening. Multiple linear regression (MLR; R2=0.9790) and support vector machine (SVM; R2=0.9114) models consistently predicted potential bioactivity of the TCM candidates against IN. The 40 ns molecular dynamics (MD) simulation showed that the TCM compounds fulfilled the drug-like criteria of forming stable complexes with IN. Atomic level investigations revealed that 9-hydroxy-(10E)-octadecenoic acid binds at an important residue A:Lys173, and Beauveriolide I formed stable interactions with the core LEDGF binding site and with Asn256 of the IN binding site on LEDGF. The TCM candidates also initiate loss of α-helices that could affect functionality of IN. Taken together, the ability of 9-hydroxy-(10E)-octadecenoic acid and Beauveriolide I to (1) form stable interactions affecting IN-LEDGF binding, and (2) have predicted bioactivity against IN suggests that the TCM candidates might be potential starting structures for developing compounds that may disrupt IN-LEDGF binding.
Keywords: human immunodeficiency virus integrase (HIV-IN), len epithelium derived growth factor (LEDGF), traditional Chinese medicine (TCM), docking, molecular dynamics (MD) simulation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
Introduction
Human immunodeficiency virus (HIV) is a type of reverse-transcription virus that utilizes host cells for survival. Following entry into the host cell, HIV will reverse transcribe its RNA into DNA, which is then integrated with host DNA at the nucleus through the action of HIV integrase (IN) . Len epithelium derived growth factor (LEDGF)/p75 is the primary binding partner of IN in host cell nuclei and plays critical roles in protecting IN from degradation and facilitating IN trafficking in
chromosomes and nucleus .
LEDGF protein is formed by 530 amino acids, and is found ubiquitously in human cells . There are two binding sites in LEDGF: the N-terminal human DNA binding site (residues 1-347) and the C terminal HIV IN binding site located from residues 348-530 . IN are proteins composed of three structural domains: the N- terminal domain (NTC), the catalytic core domain (CCD), and the C-terminal domain (CTD). NTC (residues 1-49) has four conserved amino acids that interact with zinc ion to stabilize binding of IN with viral DNA. CCD from residues 50-212 contains important active residues Asp64, Asp116, and Glu152 (also known as the DDE motif). The C-terminal domain (residues 213-288) is responsible for the non-specific binding of the IN dimer with host DNA . CCD binds with the integrase binding domain (IBD; residues 348-530) of LEDGF . Since proteosomal destruction of IN 1
2 3 4 5 6 7
8 9 10 11 12 13 14 15 16 17 18 19
occurs when IN is not bound to LEDGF, it may be possible to inhibit viral replication by blocking IN binding with LEDGF through manipulation of either CCD, IBD or
both binding sites.
The benzene derivative 4-[(5-bromo-4-{[2,4-dioxo-3-(2-oxo-2-phenylethyl)- 1,3- thiazolidin-5-ylidene]methyl}-2-ethoxyphenoxy)methyl]benzoic acid (D77) has been demonstrated to exhibit high affinity to CCD and inhibit reverse transcription of HIV- 1. D77 was reported to exert inhibitory activity by disrupting CCD-IBD interaction through competitive binding. Along this concept, we aimed to screen traditional Chinese medicine (TCM) to search for potential drug candidates from natural resources. TCM are medicinal materials used for thousands of years in Asian countries, therefore use of these TCM resources may provide safer alternatives to
synthetic drugs that often cause adverse reactions in human bodies.
Materials and Method
Protein Structure Information
The sequence of the IN was adopted from Swiss-Prot (HIV integrase_HUMAN, P12497), and the protein structure was downloaded from Protein Data Bank (PDB:
2B4J) . The IN-LEDGF complex was adopted in order to identify TCM compounds that could potentially alter the IN-LEDGF complex and elicit changes that might
induce the separation of IN and LEDGF.
Docking 1
2
3 4 5 6 7 8 9 10 11
12 13 14 15 16 17 18 19 20
TCM ligands used for docking were downloaded from TCM Database@Taiwan (http://tcm.cmu.edu.tw/). The Prepare Ligands module in Discovery Studio 2.5 (DS 2.5) was applied to prepare ligands and screen through Lipinski’s Rule of Five. The 9,029 prepared compounds that passed Lipinski’s Rule of Five were docked using the LigandFit module and the forcefield of HARvard Macromolecular Mechanics
(CHARMm) . Reported binding residues of Thr174 on A chain, and Gln95, Thr125, and Trp131 on B chain) were used to define the binding site boundaries, and the Calculate Binding Energies protocol was used to calculate binding energies. Ligands
were ranked using Binding Energy and Dock Score as the primary and secondary criteria, respectively. D77 was used as the control. Docking poses of each candidate
and D77 were evaluated with Ligplot to analyze CCD-ligand interactions.
Construction of 2D-Quantitative Structural Activity Relationship Models
Ligands with demonstrated inhibition activity against IN were adopted from Tang et al. to construct QSAR prediction models for IN inhibition (Table 1).
Descriptors for the 23 ligands and TCM candidates were calculated by DS 2.5 Calculate Molecular Property module. Descriptors in which TCM candidate values
fell within the maximum and minimum values generated for the 23 ligands were selected. Six representative descriptors were selected through Genetic Approximation (GA) algorithm calculations, and used to determine the optimum pIC50 prediction 1
2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19
model. The 23 adopted ligands were randomly divided into a training set and test set of 18 and 5 ligands, respectively. A linear Multiple Linear Regression (MLR) model incorporating the six representative descriptors was constructed using Matlab in the
following form:
Y = b0 + b1*X1 + b2*X2 + ... + b6*X6
where Y is the predicted pIC50 value, b is a constant, and X represents the descriptor.
Identical representative descriptors were also applied to construct a Support Vector Machine (SVM) model. SVM can be viewed as a multi-dimensional transformation, and can be used in non-linear regression predictions. Descriptors were first scaled to a range between -1 and 1. LibSVM was used to construct the SVM model and validity of the model was verified by cross validations inherent in the
algorithms.
Both MLR and SVM models were subjected to external validation using the randomly selected test set ligands. Validated models were then used to predict pIC50
values of the TCM ligands and D77 against IN.
Molecular Dynamics Simulation
To test protein-ligand stability, TCM candidates were docked to CCD and the complex minimized by the Steepest Descent and Conjugate Gradient for a maximum of 500 steps each using the Standard Dynamics Cascade and Dynamics (Production) package. The system was heated to the target temperature of 310K in 50 ps and 1
2 3
4 5 6 7 8 9 10 11
12 13 14
15 16 17 18 19 20
equilibrated for 200 ps. MD production time was 40 ns under constant temperature (NVT) using Berendsen weak coupling method. The temperature coupling decay time was set to 0.4 ps. MD simulation trajectories were analyzed for changes in protein and ligand root mean square deviation (RMSD), total energy, hydrogen bond (H-bond) occupancy and distance, and torsion. MD conformations were subjected to clustering to determine the representative structures of each protein-ligand complex following
MD stabilization.
Results and Discussion
Docking
Table 2 lists the Binding Energy and Dock Score for D77 and the top two TCM candidates from virtual screening. D77 (Figure 1A) has a binding energy of -85.1891 kcal/mol. Docking pose snapshots (Figure 2A) reveals that D77 forms H-bonds with A:Lys173 and B:Glu96. 9-Hydroxy-(10E)-octadecenoic acid (Figure 1B) had a predicted binding energy of -135.383kcal/mol, and formed a single H-bond with A:Lys173 (Figure 2B). Beauveriolide I (Figure 1C) was predicted to have a binding energy of -97.761kcal/mol and formed H-bonds with B:Gln94 and B:Gln95 (Figure 2C). Ligplot analysis indicated that D77 formed hydrophobic interactions with B:Gln95, and Asn367, Leu368, and Asp369 in the IBD of LEDGF (Figure 3A).
Amino acids that formed hydrophobic interactions with 9-hydroxy-(10E)- 1
2 3 4 5 6
7 8 9 10 11 12 13 14 15 16 17 18 19
octadecenoic acid (Figure 3B) were located on opposite sides of the hydrophobic tail region, forming a hydrophobic pocket in which 9-hydroxy-(10E)-octadecenoic acid was enclosed. For Beauveriolide I, hydrophobic interactions were detected not only between CCD dimers (B:Thr93, B:Thr125, and A:His171) but with IBD residues C:Asp366 and C:Asp369 as well (Figure 3C). Comparisons of our docking poses indicate that the ligands were located within the CCD reported in the literature. H- bond formation of our test candidates with B:Gln95 was consistent with previous
reports .
2D-QSAR models and bioactivity predictions
Molecular descriptors most related to pIC50 values were identified using GA.
Results of different models and their corresponding r2 values are listed in Table 3. The top model with the highest r2 of 0.9064 (implicating the highest correlation with bioactivity) is:
Zlength Shadow
Y PMI JursTASA
Kappa Zagreb
ility So
Molecular E
_ 72645
. 0 _ 0016825
. 0 0058888
. 0
1 _ 5239
. 1 31035
. 0 lub
_ 658
. 1 2568 . 3 1
Descriptions of the selected descriptors are listed in Table 4.
Figure 4 illustrates correlation of experimental and predicted pIC50 values determined using 2D-QSAR models constructed with the representative descriptors.
The high R2 values for both MLR (R2= 0.9790) and SVM (R2=0.9114) models indicate good correlation between reported experimental values and model 1
2 3 4 5 6 7
8 9 10 11 12 13 14
15 16 17 18 19
predictions. Bioactivities predicted by the MLR and SVM models for the TCM candidates and D77 are listed in Table 2. MLR predictions suggested that 9-hydroxy- (10E)-octadecenoic acid may show higher inhibitory actions against IN than D77.
Beauveriolide I was predicted to have a lower bioactivity. Alternatively, all three test ligands were predicted to have similar bioactivities by the SVM model. Regardless of prediction differences, both models suggest that TCM candidates might exhibit
bioactivity against IN.
Molecular Dynamics Simulation
Docking results are generated with a static protein and may not appropriately reflect under actual poses under physiological conditions. To account for these possible differences, we then conducted MD simulation to analyze stability of the protein-ligand complex under dynamic conditions and investigate important stabilizing interactions. Trajectories shown in Figure 5 indicate that D77 and Beauveriolide_I exhibited similar complex and ligand RMSD trajectories that were relatively stable. 9-Hydroxy-(10E)-octadecenoic acid, on the other hand, exhibited high variation which might be related to the flexibility of its hydrocarbon tail. Total energy trajectories indicated that the test ligands stabilized after 30 ns and ranged between -1,800 and -2,000 kcal/mol (Figure 5). RMSD and total energy values suggest that the complexes are stable.
1 2 3 4 5 6
7 8 9 10 11 12 13 14 15 16 17 18 19
Stability of various bonds under dynamic conditions can be examined from H- bond distance trajectories. D77 formed stable H-bonds with A:Lys173 and B:Thr124 (Figure 6A). Repetitive distance fluctuations observed between D77 and A:Lys173 were most likely due to rotation of H atoms on the lysine NH3 moiety. The H-bond with B:Thr124 stabilized after 3.36 ns. A change in benzene direction (Figure7A) might have contributed to this additional H-bond. Intriguingly, a rapid decrease in distance was observed for 9-hydroxy-(10E)-octadecenoic acid prior to the formation of stable H-bonds with A:Lys71, A:Leu172, and A:Lys173 (Figure 6B). 9-Hydroxy- (10E)-octadecenoic acid is composed of single bonds and has a highly flexible hydrocarbon tail. As shown in Figure 7B, 9-hydroxy-(10E)-octadecenoic acid undergoes large positional shifts and rotations during MD. 9-Hydroxy-(10E)- octadecenoic acid could rotate continuously until it reaches a specific conformation that enables H-bond formation and stabilizes. Beauveriolide I formed stable H-bonds with B chain residues B:Gly94, B:Gln95, and B:Thr125 at distances of approximately 2Å (Figure 6C). Contrary to D77 and 9-hydroxy-(10E)-octadecenoic acid,
Beauveriolide I did not form H-bonds with A:Lys173.
A better understanding of H-bond stability can be provided by torsion analysis (Figure 8). Torsions 1, 2, and 3 are related to H-bond formation with Thr124 by D77.
Cumulative fluctuations at torsion 1, 2, and 3 may have contributed to benzene 1
2 3 4 5 6 7 8 9 10 11 12 13 14 15
16 17 18 19
rotation enabling H-bond formation Thr124 (Figure 7A). The carboxyl group of D77 was highly dynamic (4, Figure 8A), providing a possible explanation to the H-bond distance variations recorded in Figure 6A. Large torsion fluctuations were observed for 9-hydroxy-(10E)-octadecenoic acid (Figure 8B) which are expected due to the single bond nature of the hydrocarbon tail region. Torsion angles in Beauveriolide I are smaller than 90 degrees, implying small rotational changes in H-bonds (Figure
8C).
Intriguing results can be obtained when data from docking and MD simulation are compared (Table 5). D77 and 9-hydroxy-(10E)-octadecenoic acid maintained H- bonds with A:Lys173. Beauveriolide I, on the other hand, maintained interactions with B:Gly94, B:Gln95, B:Thr125. The ability of 9-hydroxy-(10E)-octadecenoic acid to maintain interactions with A:Lys173 has important drug implications. Cherepanov et al. reported that IN mutations at A:Lys173, as well as Val165, Arg166, Leu172, Gln168, would lead to a loss of binding ability to LEDGF despite the lack of direct contact . This loss could be the result of an altered IN structure in which LEDGF no longer recognizes. Binding of 9-hydroxy-(10E)-octadecenoic acid at A:Lys173 could corrupt recognition of IN by LEDGF, hence limit successful binding. Beauveriolide I formed interactions within the core LEDGF binding region on IN (B:94-107, B:124- 133, and A: 166-178) previously reported . In particular, Beauveriolide I formed H- 1
2 3 4 5 6
7 8 9 10 11 12 13 14 15 16 17 18 19
bonds with B:Gln95 and B:Thr125. In addition, a H-bond was detected between Beauveriolide I and C:Asn367 of LEDGF. The ability of Beauveriolide I to form direct interactions at locations adjacent to direct binding residues of CCD-IBD suggest possible disruption of interface interactions A:Glu170/A:His171 with C:Asp and the hydrophobic contact of B:Leu102, B:Ala-128, B:Ala129, B:Trp132,
A:Thr174, B:Met178 with C:Ile365.
Analysis of protein-ligand complex structures help visualize structural changes elicited by the test compounds during MD, and may provide insights to possible mechanisms. In D77, α-helices in A:Gly197-A:Val201, B: Ala195-B:Thr206, and B:
Pro145-B:Asp167 were lost in the representative MD structure following stabilization (Figure 9A). Similarly, relaxation of the helix structures were also observed at A:Lys186-A:Ile208, B:Gly140-Asp167, and B:Gly197-B:Ile208 in 9-hydroxy-(10E)- octadecenoic acid (Figure 9B). α-Helices at A:Asn184-A:Gly197, B:Thr93-Glu96, and B:Phe:139-B:Gly149 were also lost when Beauveriolide I was docked into CCD (Figure 9C). Since α-helices are important structures for stabilizing protein structures, loss of these stabilizing conformations suggests changes to the CCD following binding of the TCM compounds that may alter not only its structure, but its
functionality as well.
Conclusion
1 2 3 4 5
6 7 8 9 10 11 12 13 14 15 16 17 18 19
TCM compounds 9-hydroxy-(10E)-octadecenoic acid and Beauveriolide I were predicted to have higher binding affinities and Dock Scores than the known IN agonist D77. QSAR model predicted that the TCM candidates had bioactivity against IN. Stability of TCM candidate-IN complexes further suggest that the TCM candidates have desirable drug-like characteristics. Protein-ligand interactions following stabilization show that 9-hydroxy-(10E)-octadecenoic acid binds at the important residue A:Lys173 and Beauveriolide I forms interactions with key CCD residues and Asn367 of IBD. Both TCM compounds elicit structural changes to CCD that result in the loss of important α-helices. Based on these observations, we speculate that the TCM candidates might alter interactions between IN and LEDGF, and suggest that 9-hydroxy-(10E)-octadecenoic acid and Beauveriolide I might be potential starting structures for developing compounds that may have implications in
HIV by disrupting IN-LEDGF binding.
Acknowledgements
The research was supported by grants from the National Science Council of Taiwan (NSC101-2325-B-039-001), Asia University (100-asia-56, asia100-cmu-2), and China Medical University (DMR-101-094). This study is also supported in part by Taiwan Department of Health Clinical Trial and Research Center of Excellence (DOH101-TD-B-111-004) and Taiwan Department of Health Cancer Research Center 1
2 3 4 5 6 7 8 9 10 11 12
13 14 15 16 17 18 19
of Excellence (DOH101-TD-C-111-005). We are grateful to Asia University and the
National Center of High-performance Computing for computer time and facilities.
References
Discovery Studio Modeling Environment v.2.5 (Accelerys Software Inc., San
Diego, 2009).
Bonchev, D. in Chemometrics Series Vol. 5 (ed Bawden, D.D.) (Research Studies Press Ltd., 1983).
Brooks, B.R., Brooks, C.L., Mackerell, A.D., Nilsson, L., Petrella, R.J., Roux, B., ... Karplus, M. (2009). CHARMM: The biomolecular simulation program.
Journal of Computational Chemistry, 30, 1545-1614
Chang, C.C., & Lin, C.J. (2011). LIBSVM : a library for support vector machines. ACM Transactions on Intelligent Systems and Technology, 2,
27:21-27:27
Chen, C.Y.C. (2011). TCM Database@Taiwan: the world's largest traditional
Chinese medicine database for drug screening in silico. PLoS One, 6, e15939 Cherepanov, P., Devroe, E., Silver, P.A., & Engelman, A. (2004).
Identification of an evolutionarily conserved domain in human lens epithelium-derived growth factor/transcriptional co-activator p75 (LEDGF/p75) that binds HIV-1 integrase. Journal of Biological Chemistry, 279, 48883-48892
1
2 3 4
5 6 7
8 9
10 11 12
13 14
15 16 17 18 19 20
Cherepanov, P., Ambrosio, A.L.B., Rahman, S., Ellenberger, T., & Engelman, A. (2005). Structural basis for the recognition between HIV-1integrase and transcriptional coactivator p75. Proceedings of the National Academy of Sciences of the United States of America, 102, 17308-17313
Christ, F., Voet, A., Marchand, A., Nicolet, S., Desimmie, B.A., Marchand, D., ... Debyser, Z. (2010). Rational design of small-molecule inhibitors of the LEDGF/p75-integrase interaction and HIV replication. Nature Chemical Biology, 6, 442-448
Ciuffi, A., & Bushman, F.D. (2006). Retroviral DNA integration: HIV and the
role of LEDGF/p75. Trends in Genetics, 22, 388-395
Du, L., Chen, J., Yang, L.M., Zheng, Y.T., Tang, Y., Shen, X., & Jiang, H.L.
(2008). D77, one benzoic acid derivative, functions as a novel anti-HIV-1 inhibitor targeting the interaction between integrase and cellular LEDGF/p75.
Biochemical and Biophysical Research Communications, 375, 139-144
Hall, L.H., & Kier, L.B. (1991). The molecular connectivity chi indexes and kappa shape indexes in structure-property modeling. Reviews in Computational Chemistry, 2, 367-422
Kier, L.B. (1987). Indexes of molecular shape from chemical graphs.
Medicinal Research Reviews, 7, 417-440 1
2 3
4 5 6 7
8 9
10 11 12 13
14 15 16
17 18 19
Llano, M., Delgado, S., Vanegas, M., & Poeschla, E.M. (2004a). Lens epithelium-derived growth factor/p75 prevents proteasomal degradation of
HIV-1 integrase. Journal of Biological Chemistry, 279, 55570-55577
Llano, M., Vanegas, M., Fregoso, O., Saenz, D., Chung, S., Peretz, M., &
Poeschla, E.M. (2004b). LEDGF/p75 determines cellular trafficking of diverse lentiviral but not murine oncoretroviral integrase proteins and is a component of functional lentiviral preintegration complexes. Journal of Virology, 78,
9524-9537
Maertens, G., Cherepanov, P., Pluymers, W., Busschots, K., De Clercq, E., Debyser, Z., & Engelborghs, Y. (2003). LEDGF/p75 is essential for nuclear and chromosomal targeting of HIV-1 integrase in human cells. Journal of Biological Chemistry, 35, 33528-33539
Tang, J., Maddali, K., Metifiot, M., Sham, Y.Y., Vince, R., Pommier, Y., &
Wang, Z. (2011). 3-Hydroxypyrimidine-2,4-diones as an inhibitor scaffold of
HIV integrase. Journal of Medicinal Chemistry, 54, 2282-2292
Wallace, A.C., Laskowski, R.A., & Thornton, J.M. (1995). LIGPLOT: a program to generate schematic diagrams of protein-ligand interactions.
Protein Engineering, 8, 127-134 1
2
3 4 5 6 7
8 9 10 11
12 13 14
15 16 17 18
Yen, H., Wang, G., Wang, J., & Chu, Y. (2004). Research developments in the components and pharmaceutical effects of Bombyx mori L. infected with
Beauveria bassiana (article in Chinese). SIlkworm Industry in China, 25, 86-88 Yokoi, H., Mizukami, H., Nagatsu, A., Ohno, T., Tanabe, H., & Inoue, M.
(2009). Peroxisome proliferator-activated receptor gamma ligands isolated from adlay seed (Coix lacryma-jobi L. var. ma-yuen STAPF.). Biological &
pharmaceutical bulletin, 32, 735-740
Zheng, Y., Ao, Z., Jayappa, K.D., & Yao, X. (2010). Characterization of the HIV-1 integrase chromatin- and LEDGF/p75-binding abilities by mutagenic analysis within the catalytic core domain of integrase. Virology Journal, 7, 68 1
2
3 4 5 6
7 8 9 10
11 12
1
